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ABSTRACT 


Numerical  Calculation  of  Aspiration  Efficiency  of  Aerosols 
into  Thin  Walled  Sampling  Inlets. 


Efficiency  of  aspiration  of  particles  from  a  flowing 
airstream  into  a  thin  walled  sampling  inlet  is  accurately 
predicted  using  a  numerical  model.  The  model  combines  the 
Boundary  Integral  Equation  Method  for  predicting  the  air 
velocity  field  into  the  inlet  with  an  analytical  solution  to 
the  particle  equations  of  motion  to  describe  particle 
trajectories.  A  two  and  a  three  dimensional  model  are 
examined.  Results  for  the  three  dimensional  model  correlate 
well  with  an  empirical  model  of  aspiration  efficiency.  The 
method  can  be  generalized  for  a  wide  range  of  airstream, 
sampling  inlet,  and  particle  conditions. 


Author:  Kevin  M.  Boyle,  Capt,  USAF 

Year : 1991 

No  of  Pages:  119 

Degree:  Master  of  Science  in  Environmental  Engineering 
Name  of  Inst:  University  of  North  Carolina  at  Chapel  Hill 


91-07319 


TABLE  OF  CONTENTS 


I  .  BACKGROUND . 1 

1 1  .  THEORY . 7 

III.  METHOD . 13 

IV.  RESULTS . 31 

V.  DISCUSSION . 39 

VI.  CONCLUSION . 48 

VII.  FUTURE  WORK . 49 

VIII.  APPENDIX . 53 

IX.  REFERENCES . 118 


iii 


ACKNOWLEDGEMENTS 


I  would  like  to  thank  Dr.  Michael  R.  Flynn,  Department 
of  Environmental  Sciences  and  Engineering,  for  his  guidance 
and  constant  encouragement  throughout  the  course  of  this 
project.  I  am  also  grateful  to  Dr.  David  Leith  and  Dr. 
Russell  W.  Wiener  for  their  comments  and  suggestions 
regarding  the  content  and  preparation  of  this  report. 

I  extend  my  sincerest  appreciation  to  my  wife,  Ellen, 
and  daughter,  Sarah,  for  their  understanding  and  support 
throughout  this  project. 

I  would  also  like  to  thank  the  United  States  Air  Force 
for  providing  the  opportunity  and  financial  support  needed 
to  complete  my  program  of  study. 

This  research  was  funded  by  the  U.S.  Environmental 
Protection  Agency  under  grant  NO.  44213. 


NUMERICAL  CALCULATION  OF  INERTIAL  ASPIRATION  EFFICIENCY  OF 


AEROSOLS  INTO  THIN-WALLED  SAMPLING  INLETS 

by 

Kevin  M.  Boyle 

Unbiased  sampling  of  airborne  particulate  from  a 
flowing  stream  requires  that  the  size  distribution  and 
concentration  of  aerosol  collected  be  identical  to  that  of 
the  aerosol  in  the  free  stream.  Sampling  errors  occur 
during  aspiration  of  the  aerosol  from  the  free  stream  to  the 
face  of  the  inlet  and  during  transmission  of  the  aerosol 
along  the  sample  tube.  Additional  losses  or  gains  may  occur 
due  to  particle  bounce  from  the  front  edge  of  the  sample 
tube.  In  this  paper,  a  numerical  model  for  determining  the 
aspiration  component  of  overall  sampling  efficiency  for  any 
arbitrarily  shaped  thin  edged  inlet  is  presented. 

BACKGROUND 


Sampling  errors  due  to  aspiration  occur  when  particles 
fail  to  enter  the  sampling  inlet  with  the  air.  Particles 
have  inertia  which  may  prevent  them  from  accelerating, 
decelerating,  or  changing  their  direction  fast  enough  to 
move  with  the  air  (Agarwal  and  Liu,  1980).  Particles  also 
possess  additional  motion  due  to  gravity.  Aerosol 
aspiration  efficiency,  Ea,  is  defined  in  equation  (1)  as  the 
ratio  of  particle  concentration  entering  the  sampling  inlet. 
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Ci,  to  the  particle  concentration  at  an  undisturbed  upstream 
location,  Co,  (Belyaev  and  Levin,  1972). 

Ea  =  Ci  (1) 

Co 

One  effect  on  aspiration  efficiency  would  be  loss  or 
gains  of  particles  due  to  bounce  off  the  front  of  the 
sampling  inlet.  By  modeling  the  system  as  a  thin  walled 
sampling  inlet,  particle  rebound  effects  are  assumed 

negligible  and  do  not  contribute  to  the  overall  aspiration 
efficiency  (Belyaev  and  Levin,  1974). 

Both  numerical  modeling  and  physical  experiments  have 
attempted  to  describe  aspiration  efficiency.  Belyaev  and 
Levin  (1972,  1974)  used  flash  illumination  to  photograph 
particle  trajectories  and  develop  equations  to  predict 
inertial  aspiration  efficiency.  Durham  and  Lundgren  (1980), 
developed  aspiration  equations  from  experimental  data  for 
aspiration  efficiency  during  anisokinetic  sampling  with 
nozzle  misalingment  and  sampler  inlet  velocity  different 
from  oncoming  wind  velocity.  Their  equations  account  for 
the  projected  area  of  the  nozzle  when  it  is  misaligned. 
Jaysekera  and  Davies  (1980)  and  Davies  and  Subari  (1982) 
also  performed  experiments  and  derived  aspiration  equations 
which  were  similar  to  the  equations  of  Belyaev  and  Levin. 

Extensive  measurements  of  sampling  efficiency  were  made 
by  Tufto  and  Willeke  (1982),  and  Okazaki,  Wiener,  and 
Willeke  ( 1987a, 1987b, 1987c )  for  both  isoaxial  and  non- 
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isoaxial  aerosol  sampling  conditions.  Based  on  this  data 
and  work  cited  above,  Hangal  and  Willeke  (1990)  proposed  a 
unified  model  to  predict  aspiration  efficiency,  Ea.  The 
general  expression  ,  equation(2),  is: 

Ea  =  1  +  (R  cos  e  -  l)f  (2) 

where  R  is  the  velocity  ratio,  6  is  the  sampling  angle  and  f 
is  an  inertial  parameter.  For  the  case  of  0°  sampling  angle 
the  inertial  parameter,  f,  is  based  on  the  work  of  Belyaev 
and  Levin  (1974)  and  is  shown  in  equation  (3). 

f  =  1  -  1  /[I  +  (2  +  0.617/R)  Stk]  (3) 

R  =  Uo/Ui 

Uo  is  the  free  stream  velocity 
Ui  is  the  inlet  velocity 
Stk  is  the  Stokes  number 

For  0°  to  60°  sampling  angle  the  parameter  is  based 
on  the  work  of  Durham  and  Lundgren  (1980),  equation  (4). 

f  =  [1  -1/[1  +(2  +  0.617/R)stk' ] ] 

[1  -  1/[1  +  O.SSstk'  exp( 0 . 25stk ' ) ] ] / 

[1  -  1/[1  +  (2.617)stk’ ] ]  (4) 

where 

stk'  =  Stk  exp  (0.0226)  (5) 

To  model  the  movement  of  particles  in  air  it  is 
necessary  to  model  the  airflow  field  and  derive  a  method  to 
track  particles  placed  in  the  air  stream.  (Agarwal  and  Liu, 
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1980).  Numerical  models  of  air  flow  fields  and  particle 
trajectory  paths  have  been  developed  by  Agarwal  and  Liu 
(1980),  Addlesee  (1980),  and  Dunnett  and  Ingham  (1986, 

1988).  Agarwal  and  Liu  determined  the  theoretical  flow  field 
into  a  thin  walled  sampling  inlet  considering  both  particle 
inertia  and  gravitational  settling.  A  numerical  calculation 
procedure  was  used  to  solve  the  axisymmetric  Navier-Stokes 
equations  and  particle  motion  equations. 

Addlesee,  (1980)  studied  the  intake  efficiency  of  a 
Casella  cascade  impactor.  The  impactor  was  modeled  as  a  two 
dimensional  infinite  slot  with  no  end  effects.  A  solution 
to  the  flow  field  was  developed  by  assuming  that  viscous 
forces  were  not  significant  and  calculating  stream  lines 
using  potential  theory.  From  the  flow  field  developed,  the 
particle  equations  of  motion  were  solved  assuming  that 
Stokes  law  is  applicable  (Re  <  1).  The  equations  were  put 
into  finite  difference  form  and  solved  by  relaxation. 

Dunnett  and  Ingham  (1986)  developed  a  mathematical 
theory  to  predict  aspiration  efficiency  of  two  dimensional 
blunt  bodies  using  the  Boundary  Integral  Equation  Method. 
Once  the  velocity  field  of  the  system  was  known,  the 
particle  equations  of  motion  are  solved  assuming  that  the 
Reynolds  number  is  always  less  than  1.0  and  thus  Stokes' 
drag  could  be  applied.  The  equations  were  put  in  finite 
difference  form  and  solved  numerically.  In  1988  they 
expanded  the  two  dimensional  solution  to  an  axisymmetric 
blunt  body  sampler  and  numerically  predicted  aspiration 


efficiency  by  particle  tracking  using  finite  difference 
techniques . 
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Flynn  and  Miller  (1989)  used  the  boundary  integral 
equation  method  to  model  air  flow  into  local  exhaust  hoods. 
Their  solutions  follow  the  same  logic  as  Dunnett  and  Ingham 
and  give  good  agreement  with  analytic  solutions  for  flow 
into  an  infinitely  flanged  rectangular  hood  with  a  constant 
velocity  at  the  hood  face,  and  with  other  empirical  models. 
The  advantage  of  their  model  was  that  it  reduced  a  3 
dimensional  air  flow  solution  to  a  2D  area  solution  along 
the  boundary  of  a  specified  domain. 

This  research  employs  the  model  of  Flynn  and  Miller  to 
estimate  the  air  flow  into  a  sampler  inlet,  and  then  applies 
the  solution  for  the  particle  equations  of  motion  developed 
by  Alenius  (1989)  to  determine  aspiration  efficiency.  The 
particle  motion  equations  are  described  in  the  THEORY 
section  which  follows. 

The  referenced  numerical  models  to  calculate  aspiration 
efficiency  were  limited  by  the  geometry  of  the  solution 
domain  and  employed  finite  difference  techniques  to  solve 
the  particle  equations  of  motion.  By  combining  the  3D 
airflow  solution  of  Flynn  and  Miller  with  the  analytic 
solution  to  the  equations  of  motion  proposed  by  Alenius 
(1989),  this  model  can  calculate  particle  trajectories  into 
a  thin  edged  inlet.  This  model  has  a  significant  advantage 
over  other  models  as  it  does  not  require  finite  difference 
discretization  to  calculate  particle  trajectories  and  has 


the  flexibility  to  be  developed  for  any  arbitrarily  shaped 
domain  and  sampling  inlet  over  a  wide  range  of  particle  size 
and  velocity  conditions. 
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THEORY 

To  model  particle  flow,  a  method  for  determining  the 
characteristics  of  the  flow  field  into  an  arbitrarily  shaped 
and  oriented  sampling  inlet  is  required.  Once  the  flow 
field  is  described,  particles  can  be  placed  into  the  field, 
assigned  initial  conditions,  and  tracked.  From  the 
trajectories,  an  estimate  of  the  aspiration  efficiency  can 
be  made  and  compared  with  previous  models  or  experimental 
data . 


The  boundary  integral  equation  method  (BIEM)  is  a 
method  of  solving  a  potential  flow  problem  for  any 
arbitrarily  shaped  domain.  The  assumptions  in  the 
development  are  inviscid,  incompressible,  irrotational  flow 
which  leads  to  the  following  expressions: 

grad-V  =  0  (continuity  equation  for  (6) 
incompressible  flow  ) , 

and ; 

curl  V  =  0  (irrotational  flow).  (7) 
Irrotational  vector  fields  possess  a  potential  function  such 
that  the  vector  is  equal  to  the  gradient  of  the  potential , 
d,  and  thus 


V  =  grad  d . 

By  combining  the  above  equations 


(8) 


grad-V  =  grad -(grad  d)  or 
grad  ^0  =  0 


(9) 
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Equation  (9)  is  Laplace's  equation.  By  applying  the 
divergence  theorem  and  Green's  second  identity,  an  equation 
can  be  derived  which  relates  the  velocity  potential  an  any 
singular  point  (P)  to  the  potential  and  normal  derivatives 
of  potential  along  the  boundary  of  a  specified  domain. 

The  derivation  (Liggett  and  Liu,  1983)  yields  the 
following  equations  which  form  the  basis  of  the  Boundary 
Integral  Equation  Method  for  two  and  three  dimensional 
domains . 


2D: 

ad(P)  = 

3D: 

-ad(P)  = 

Sr, 


[((d/R)  6R/6n) 
[(d(6(l/R)/6n) ) 


((InR)  6d/6n)]  dS 
( ( l/R)6d/6n) ]  dA 


(10) 

(11) 


6 

66/6X1 

a 

R 

r 

dS 

dA 


velocity  potential 

velocity  normal  to  the  boundary 

solid  angle  about  (P) 

distance  from  singular  point  (P)  to 

a  point  on  the  boundary. 

closed  curve  surrounding  and  defining 

the  2D  domain 

closed  surface  surrounding  and  defining 
the  3D  domain 
differential  element  on  P 
differential  area  on  Q. 


The  above  equations  can  be  written  in  discrete  form  in 
either  two  (Liggett  and  Liu,  1983)  or  three  dimensions 
(Flynn  and  Miller,  1989).  After  specifying  either  6  or 
66/6n  at  each  node  on  a  discretized  boundary  of  N  nodes, 
equation  (10)  or  (11)  is  written  for  each  node  by 
successively  placing  the  singular  point  (P)  at  each  node  on 
the  boundary.  This  results  in  a  system  of  N  equations  in  N 
unknowns  which  are  solved  simultaneously  to  give  the  unknown 
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value  of  0  or  60/6n  at  each  node  on  the  boundary.  Once  the 
boundary  values  are  known,  (P)  can  be  located  at  any 
internal  location  and  the  left  hand  side  of  equation  (10)  or 
(11)  differentiated  to  give  the  velocity  at  the  internal 
point . 

Once  the  velocity  field  is  determined,  particles  can 
be  placed  in  the  air  stream  and  their  motion  tracked.  This 
model  uses  an  analytical  solution  to  the  particle  equations 
of  motion  developed  by  Alenius  (1989).  The  general  form  of 
the  equation  is: 

d2r/dt2  -  [Q(Re)/T(u-(dr/dt))]-g  =  0  (12) 

where  bold  face  variables  denote  vector  quantities. 

In  the  above  equation,  the  following  notation  is  used: 
d2r/dt2  (=dv/dt)  is  the  acceleration  of  the  particle 

Q(Re)  is  a  dimensionless  factor  based  on  the 
Reynolds  number  and  the  drag 
coefficient 

T  is  the  particle  relaxation  time  in  air 

u  is  the  air  velocity  at  the  particle 
1 ocation 

dr/dt=v  is  the  velocity  of  the  particle 

g  is  the  earth's  gravity 

0  is  the  null  vector 

r  is  the  particle  position  vector 

The  general  equation  presented  is  not  analytically 
solvable.  The  reason  is  that  the  factor  Q(Re)  and  the  air 
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velocity  u  vary  along  the  trajectory  of  the  particle  and 
thus  indirectly  with  time  t.  If  the  air  velocity  at  the 
points  where  the  particle  is  successively  located  can  be 
described  as  a  linear  function  of  the  time  that  has  elapsed 
and  Q(Re)  is  assumed  to  be  constant  over  the  time  interval, 
an  analytical  solution  to  the  general  equation  can  be 
obtained.  This  condition  is  approximately  met  if  the  motion 
of  the  particle  is  successively  calculated  using 
sufficiently  small  time  intervals.  A  complete  derivation  is 
presented  in  the  appendix.  The  analytic  solution  is: 

v(t+t’)  =  Uat ’+Ub-(Ub-v(t) )exp-[Q(Re)t ' /t]  fl3) 

r(t+t')  =  r(t)+Ua/2t ' 2+Ubt ' 

-(Ub-v(t) )T/Q(Re) (l-expQ(Re)t' /T)  (14) 

where  t'  is  the  size  of  the  time  step  and  Ua  and  Ub  have 
been  added  to  simplify  the  appearance  of  the  equation: 

Ua  =  (u(t)  -u(t-t' ))/t'  (15) 

Ub  =  u(t)-(Ua-g)T/Q(Re) .  (16) 

Particles  can  now  be  tracked  in  any  domain  in  which  the 
air  field  is  known,  and  the  efficiency  of  aspiration 
determined.  Aspiration  efficiency  is  calculated  as  the 
ratio  of  particles  entering  the  inlet  to  particles  in  the 
undisturbed  flow  field  far  upstream,  equation  (1).  Figure  1 
shows  a  conceptual  stream  tube  of  particles  flowing  in  an 
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Freestream ; 


Figure  1.  Conceptual  Stream  Tube.  The  flux  of  particles 
through  the  upstream  plane  equals  that  of  particles  passing 
the  face  of  the  ssonpler  inlet. 
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airstream  of  given  velocity,  Uo,  into  an  inlet  sampling  at  a 
constant  velocity,  Ui .  The  boundaries  of  the  stream  tube 
represent  all  particle  trajectories  which  just  enter  the 
tube.  These  trajectories  are  called  critical  trajectories. 
Particles  located  outside  this  tube  would  not  be  drawn  into 
the  inlet.  The  flux  of  particles  passing  through  the 
upstream  plane  of  the  stream  tube  is  particle  concentration, 
Co,  multiplied  by  the  stream  tube  cross  sectional  area,  Ao , 
times  the  freestream  velocity,  Uo .  By 

equating  the  flux  of  particles  through  the  plane  of  the 
stream  tube  upstream  from  the  inlet  to  that  through  the 
plane  of  the  sampler  face,  Ci  Ai  Ui,  the  expression  in 
equation  (17)  is  obtained. 

Co  Ao  Uo  =  Ci  Ai  Ui  (17) 

From  equation  (1),  Ea  =  Ci/Co,  thus: 

Ea  =  Ao  Uo/Ai  Ui  (18) 

In  equation  (18),  the  characteristics  of  the  sampler 
inlet,  Ai  and  Ui,  are  known  as  well  as  the  free  stream 
velocity,  Uo .  By  determining  the  coordinates  of  the 
critical  trajectory  start  points  in  the  free  stream,  an 
estimate  of  Ao  and  Ea  can  be  made.  This  can  be  done  by 
iteratively  assigning  particle  locations  in  the  free  stream 
and  tracking  each  particle  in  the  domain  until  the  critical 
trajectories  are  located. 
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METHOD 


The  model  was  developed  using  a  simple  two  dimensional 
problem  which  was  then  expanded  to  a  more  complex  three 
dimensional  case.  Both  isoaxial  and  non-isoaxial  inlet 
orientations  were  modeled.  Non-isoaxial  conditions 
consisted  of  orienting  the  inlet  at  a  pitch  angle  of  30®  to 
the  freestream. 

Two  Dimensional  Model 

The  two  dimensional  problem  was  modeled  as  an  infinite 
slot.  The  coordinate  system  to  visualize  in  two  dimensions 
is  an  x-y  cartesian  system  with  the  origin  located  at  the 
center  of  the  inlet.  The  positive  x  direction  extends  along 
the  centerline  of  the  hood  face  and  the  y  direction  is  up 
and  down,  with  gravity  acting  in  the  -y  direction. 

The  air  flow  into  the  slot  was  modeled  using  the 
boundary  integral  equation  method  which  provided  a  solution 
everywhere  on  the  boundary  of  the  domain  for  the  potential 
and  the  normal  derivative  of  potential .  From  the  boundary 
solution,  the  exact  velocity  components  at  any  internal 
point  were  calculated. 

A  FORTRAN  computer  program  for  the  two  dimensional 
boundary  integral  equation  solution,  written  by  Liggett  and 
Liu  (1983),  was  used  for  the  2D  solution.  The  program, 
called  GM8 ,  returned  the  boundary  values  and  the  local  air 
velocity  components  within  the  domain  needed  to  calculate 
particle  trajectories.  GM8  was  combined  with  a  subroutine 
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written  to  solve  the  Alenius  solution  to  the  equations  of 
motion,  (11)  and  (12),  and  track  particles  within  the 
domain.  Figure  2  shows  the  flow  of  logic  in  the  program. 

The  boundary  of  the  domain  in  the  simple  2-D  case  is 
shown  in  figure  3.  The  origin  of  the  coordinate  axes  are 
located  at  the  center  of  the  sampling  inlet.  The  upper 
lower,  front  and  rear  boundaries  of  the  domain  are  at  a 
distance  of  20  inlet  diameters  from  the  inlet  face.  The 
remaining  boundary  consists  of  a  rectangular  indentation  to 
represent  the  sampling  nozzle  of  diameter  equal  to  its 
width . 

In  figure  3,  the  boundary  conditions  are  marked.  The 
upper,  lower,  and  rear  boundary  are  situated  far  enough  away 
from  the  sampling  inlet  to  have  potential  due  exclusively  to 
the  free  stream  velocity,  or: 

(b  =  Uox(x)  +  Uoy(y)  (19) 

where  Uox  and  Uoy  are  the  magnitude  of  the  free  stream 
velocity  components  and  x,y  are  the  location  coordinates. 

The  boundary  condition  on  the  line  opposite  the  hood  face  is 
6d\6n  =  Uox,  the  free  stream  velocity  normal  to  the  surface. 
On  the  upper  and  lower  edges  of  the  sampling  inlet  the 
boundary  conditions  are  66\6n  =  0,  no  flow  through  the  walls 
of  the  nozzle.  At  the  face  of  the  inlet  60\6n  is  equal  to 
sampling  velocity,  Ui . 

For  the  situation  where  an  inlet  pitch  angle  was 
introduced,  the  domain  was  not  rotated  30*^  and  a  new 
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Figure  2.  Flow  diagram  for  the  two  dimension  and  3 
dlmensioal  program  logic. 
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f  =  Uox(x)  -i-  Uoy(y) 


^sUoxOO 
+  Uoy(y) 


UoxOO 
+  Uoy(y) 


Flgrure  3.  2D  model  boundary  with  coordinate  axes,  boundary 
condtlons,  and  sample  discretization  shown.  0  =  potential, 
60/6n  *s  normal  derivative  of  potential  (air  velocity) . 
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boundary  calculated.  Rather  boundary  conditions  were 
adjusted  to  represent  the  freestream,  Uo ,  as 
the  resultant  of  the  free  stream  normal  to  the  front 
boundary,  Uox,  and  a  crossdraft,  Uoy.  This  required  a 
coordinate  system  conversion  where  the  gravity  vector  in  the 
X  and  y  direction  were  -g  sin(e)  and  -g  cos(©)  respectively, 
©  was  the  angle  of  misalignment  and  g  the  magnitude  of  the 
gravity  vector.  The  origin  of  this  coordinate  system  was 
still  centered  at  the  inlet  face. 

The  solution  required  the  domain  to  be  discretized  into 
nodal  elements  as  shown  in  figure  3.  Much  finer 
discretization  than  presented  in  the  example  is  required 
because  numerical  accuracy  during  calculation  of  internal 
points  suffers  when  the  internal  point  is  within  1  element 
length  of  the  boundary. 

The  initial  conditions  assigned  to  the  particles  in  the 
free  stream  are  that  they  are  falling  at  their  terminal 
settling  velocities  and  moving  in  the  air  stream  (Addlesee, 
1980)  at  the  velocities  calculated  by  the  BIEM  solution. 
Terminal  settling  velocity  (Vts)  was  calculated  using 
Stokes'  Law: 

Vts  =  T  g  (20) 

where  t  is  the  particle  relaxation  time  (Reist,  1984)  and  g 
is  the  acceleration  due  to  gravity. 

In  two  dimensions,  it  was  necessary  to  define  the 
location,  in  the  free,  stream  of  the  upper  and  lower 
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critical  trajectories  to  calculate  aspiration  efficiencies. 
This  is  because  in  two  dimensions  equation  (18)  reduces  to: 

Ea  =  do  Uo/di  Ui  (21) 

where  di  is  the  inlet  diameter,  and  do  is  the  distance 
between  the  upper  and  lower  critical  trajectories.  Internal 
air  velocity  coordinates  were  calculated  using  BIEM,  then 
the  particle  position  and  velocity  were  predicted  using  the 
analytical  solution  and  the  process  repeated.  The  critical 
trajectories  were  determined  by  selecting  a  point  far 
upstream  and  tracking  a  particle  from  this  point  either  into 
(capture)  or  past  (miss)  the  inlet.  Once  a  missed  and 
captured  trajectory  were  found  the  average  of  the  two 
starting  points  determined  the  next  point.  Convergence 
occurred  when  the  distance  between  the  particle  trajectory 
start  point  which  is  captured  and  that  which  escapes  is  less 
than  a  preset  value.  The  program  was  designed  to  locate 
both  an  upper  and  lower  critical  trajectory  corresponding  to 
the  upper  and  lower  edges  of  the  sampling  inlet.  The 
distance,  s,  between  the  upper  and  lower  critical 
trajectories,  do,  was  determined  using  the  distance  formula, 
equation  ( 22 ) . 

s  =  ((Xi-X:)2  +  (y:-yr)2  +  (z,-z  )2)^  (22) 

Three  Dimensional  Model 

In  three  dimensions,  Uo,  Ai ,  and  Ui  are  given  and  Ao  is 
determined  to  calculate  aspiration  efficiency,  equation 
(18).  Ao  is  the  area  of  the  trajectory  tube  at  an  upstream 
location  where  the  air  flow  is  not  influenced  by  the 
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sampling  velocity.  This  area  represents  the  boundary  of  the 
stream  tube  at  an  upstream  location  in  which  particles  will 
enter  the  sampling  inlet  (Fig  1). 

The  shape  of  the  three-dimensional  domain  is  a  square 
sided  volume  of  half-width  equal  to  20  inlet  widths. 
Extending  from  the  rear  wall  of  the  volume  into  the  center 
of  the  volume  for  a  distance  of  20  inlet  widths  is  a 
sampling  inlet  with  a  square  cross  section.  An  equal  area 
square  is  used  to  model  a  circular  inlet  to  simplify  the  3D 
boundary  discretization..  The  origin  of  the  cartesian 
coordinate  system  is  located  at  the  center  of  the  inlet. 

The  X  direction  extends  away  from  the  hood  face  along  its 
centerline.  The  z  coordinate  represents  the  height  of  the 
domain  and  the  y  coordinate  represents  the  width. 

For  boundary  conditions,  the  potential  is  specified  on 
the  front,  rear,  top,  bottom  and  side  faces  of  the  3D 
domain.  The  velocity  is  specified  over  the  inlet  face  as 
Ui,  and  at  the  walls  of  the  inlet  tube  as  0  cm/sec.  As  in 
2D,  the  boundaries  of  the  3D  domain  are  situated  far  enough 
away  from  the  inlet  to  have  potential  due  exclusively  to  the 
free  stream: 

0  =  Uox(x)  +  Uoz(z)  (23) 

Here,  Uox  represents  the  velocity  perpendicular  to  the  front 
face  of  the  domain  and  Uoz  is  the  crossdraft  vector  used 
when  the  sampling  angle  did  not  equal  zero.  As  in  2D,  the 
gravity  components  in  the  coordinate  system  must  be  adjusted 
when  a  crossdraft  is  used  to  model  the  sampling  angle. 
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The  flow  field  is  defined  using  the  same  theory  as  in 
two  dimensions  but  using  a  FORTRAN  program  written  by  Flynn 
and  Miller  (1989)  to  solve  the  BIEM  problem  for  three 
dimensions . 

Particle  trajectories  were  calculated  in  the  same 
manner  as  in  two  dimensions  with  BIEM  providing  the  local 
air  velocity  components  at  the  particle  location.  The 
particles  were  assigned  initial  conditions  as  described  in 
the  2D  section  and  the  equations  of  motion,  with  a  third 
dimension  added,  were  solved  to  predict  the  particle's  next 
position  and  velocity.  The  program  logic  is  the  same  as 
that  shown  in  figure  2. 

Because  the  domain  is  now  in  three  dimensions,  the 
region  cannot  be  discretized  into  lines  between  nodes, 
rather,  the  boundary  must  be  divided  into  areas.  Flynn  and 
Miller  (1989)  used  triangular  elements  on  the  boundary. 
Examples  of  the  discretization  input  files  for  the  boundary 
are  presented  in  the  appendix. 

The  critical  trajectories  in  3D  were  determined  by 
selecting  a  point  far  upstream  and  tracking  a  particle  from 
this  point  either  into  the  inlet  (capture)  or  out  of  the 
domain  (miss).  The  start  point  is  a  specified  number  of 
inlet  diameters,  XDIN,  in  front  of  the  inlet  and  in  the 
direction  of  the  resultant  velocity  vector  in  the  x-z  plane. 
For  0°,  this  point  is  along  the  x  axis  at  x=XDIN,  y=0,  z=0, 
because  Uo  =  Uox  for  this  condition.  For  all  other  angles 
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the  first  trajectory  start  point,  figure  4,  is  calculated 
as : 

X  =  XDIN  cos(e) 

y  =  0.0  (24) 

z  =  XDIN  sin(e) 

The  rationale  for  choosing  this  start  point  is  to 
obtain  a  capture  trajectory  around  which  starting  points  for 
the  critical  trajectories  can  be  centered.  If  this  first 
point  does  not  yield  a  capture,  the  start  point  is 
incrementally  adjusted  up  or  down  along  a  line  perpendicular 
to  the  free  stream  and  the  trajectory  is  recalculated.  The 
process  is  repeated  until  a  first  capture  identified. 

Once  a  capture  trajectory  is  found,  a  series  of  start 
points  at  which  to  begin  critical  trajectories  are  defined 
at  a  distance  far  enough  away  from  the  coordinates  of  the 
first  capture  to  have  a  high  probability  of  a  miss.  These 
points  are  located  a  specified  distance,  AA ,  from  the 
coordinates  of  the  first  capture  (xl,yl,zl)  in  incremental 
degree  steps,  a,  around  a  semicircle  in  the  plane 
perpendicular  to  the  free  stream,  figure  5.  Critical 
trajectories  for  only  one  half  the  area  were  calculated  due 
to  the  symmetry  of  the  problem. 

The  equations  for  these  points  are; 

X  =  xl  +  (AA  sin(e)  sin(a)) 
y  -  AA  cos  (a)  (25) 

z  =  zl  +  ( AA  cos(©)  sin(a)) 
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Figure  4.  Location  of  initial  trajectory  start  point. 
(X1,Y1,Z1)  is  initially  placed  at  a  distance  XDIN  in  front 
of  the  inlet  along  a  line  perpencular  to  Uo  in  the  X-Z 
plane . 
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■  Trajectory  start  points 

X  =  x1  +  (AA  sin(O)  sln(o)) 
y  =  AA  cos(a) 
z  =  z1  +  (AA  cos(0)  sin(a)) 
0  =  angle  of  misalignment 


Figure  5.  Trajectory  start  points  are  located  on  a 
semicircle  centered  around  (Xl.Yl.Zl)  at  degree  increments, 
a,  along  rays  a  distance  AA  from  the  center. 

(-90®  <.  a  <  90®) 
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If  the  start  points  defined  in  equation  (25)  result  in 
a  capture,  the  start  point  is  incrementally  adjusted  farther 
from  the  initial  capture  point  and  the  trajectory  is 
recalculated  until  a  miss  identified. 

Each  first  capture  and  initial  miss  coordinate 
combination  results  in  a  ray  the  end  points  of  which  define 
the  upper  and  lower  bounds  of  the  critical  trajectory  and 
all  rays  are  contained  in  the  same  plane  perpendicular  to 
the  free  stream.  From  this  point,  the  starting  coordinates 
of  the  last  miss  and  last  capture  trajectory  are  averaged  to 
determine  the  coordinates  of  the  next  start  point  on  the 
ray.  The  process  is  repeated  until  the  convergence  criteria 
was  satisfied.  Critical  trajectories  are  found  for  each 
starting  point  around  the  semicircle. 

The  output  of  the  particle  tracking  routine  is  a  series 
of  points  representing  the  critical  trajectories  which 
define  Ao .  The  area,  Ao ,  is  calculated  numerically  with  a 
FORTRAN  program  using  the  trapezoid  rule.  From  Ao ,  Ea  is 
calculated  using  equation  (18). 

The  models,  both  2D  and  3D,  were  verified  by  comparing 
the  computed  aspiration  efficiency  with  the  unified 
aspiration  model  proposed  by  Hangal  and  Willeke  (1990), 
equation  ( 2 )  . 

The  Programs 

The  two  dimensional  program  was  compiled  and  run  using 
Microsoft  Fortran  version  2.  The  entire  aspiration 
efficiency  calculation  was  run  using  one  program.  This 
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program  initially  calculates  the  boundary  solution  for  the 
domain  then  calls  subroutines  to  calculate  particle  velocity 
and  position,  and  internal  air  velocity  components  at  the 
particle's  location.  All  2D  simulations  were  run  on  a  Dell 
386  personal  computer  and  took  approximately  15  minutes  to 
complete.  Each  program  run  calculates  approximately  30 
separate  trajectories  with  approximately  100  points 
calculated  per  trajectory. 

Figure  6  shows  an  example  of  the  output  of  the  2D 
program  with  particle  positions  plotted  as  trajectories. 

The  particle  trajectories  were  plotted  using  a  simple 
routine  written  in  BASIC  (see  appendix).  The  distance 
between  the  upper  and  lower  critical  trajectories  was  the 
basis  of  determining  the  aspiration  efficiency  as  defined  in 
equation  (22).  A  commented  version  of  the  2D  program  is 
contained  in  the  appendix. 

The  3D  program  is  much  more  complex  than  the  2D  version 
because  of  the  addition  of  a  third  dimension  and  the 
resulting  complexities  of  solving  the  BIEM  equations. 

Because  of  this,  3D  simulations  were  run  on  a  Convex  Super 
Computer.  The  3D  simulations  were  run  using  two  FORTRAN 
programs.  The  first  program,  4PTGQ.F,  is  the  BIEM  boundary 
solution  program.  The  input  for  this  program  was  a  boundary 
node  and  element  file  for  each  combination  of  sampling 
velocity,  Ui ,  wind  velocity,  Uo,  and  sampling  angle,  9.  The 
input  file  contained  367  nodes,  N,  leading  to  N  linear 
equations  and  N  unknowns  which  were  solved  using  a  four 
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Figure  6.  Sample  of  2D  trajectory  plots.  Trajectories 
originate  in  the  plane  of  the  freestream  perpendicular  to 
Uo. 
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point  gaussian  quadrature.  This  program  took  5.5  minutes  to 
run  on  the  Convex  and  returned  potential  and  velocity  values 
for  each  node  on  the  boundary  of  the  domain. 

Once  the  output  to  4PTGQ.F  was  received,  it  was 

imported  as  input  to  the  3D  particle  tracking  program, 

TRK15.F.  This  program  consists  of  a  driver  to  specify 

particle  location,  test  for  convergence,  and  call 

subroutines  to  perform  the  BIEM  internal  air  velocity, 

particle  position,  and  particle  velocity  calculations. 

TRK15.F  was  designed  to  calculate  13  trajectories  at  15°  (a 

=  15°  in  equationkk  (25))  increments  around  a  semicircle 

surrounding  the  initial  capture  point  (see  figure  5).  This 

leads  to  the  calculation  of  approximately  10  trajectories 

per  starting  point,  over  the  13  starting  points  specified  by 

the  degree  increment,  a,  with  20-30  points  along  each 

trajectory  to  determine  the  fate  of  the  particle.  A  typical 

combination  of  program  variables  shown  below  took  68  minutes 

of  computer  time  on  the  Convex: 

Uo  =  500  cm/sec 

Ui  =  250  cm/sec 

particle  diameter  =  5  pm 

e  =  30  degrees 

time  step  =  0.0005  seconds 
convergence  criteria  =  0.05  cm 

The  output  of  TRK15.F  is  a  series  of  points  which 
defines  a  region  in  the  freestream  inside  of  which  all 
particles  are  aspirated  to  the  face  of  the  sampling  inlet. 
Figures  7  and  8  show  an  example  of  the  3D  model  output. 


3D  isoaxial  Sample  Output 


z 


Plot  of  critical  trajectories 
in  freestream,  view  AA. 


Figure  7.  Plot  of  critical  trajectory  start  points  in  their 
starting  plane.  The  area  inside  the  curves  is  the  estimate 
of  Ao  in  equation  (14).  The  area  is  assumed  symmetrical 
about  \  *  0.0, 
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3D  Non-isoaxiai  Sample  Output 
Z 


0  02  0.4  0.6  O.B  1  12 


Plot  of  critical  trajectories 
in  freestream,  view  AA. 


Figure  8.  Plot  of  critical  trajectory  start  point  in  their 
starting  plane.  The  area  inside  the  curves  is  the  estimate 
of  Ao  in  equation  (14).  Z*  ,  Y'  represent  the  local 
coordinate  system  centered  at  (Xl.Yl.Zl). 


30 


Figure  7  is  an  example  of  output  for  an  isoaxial  set  of 
conditions.  The  plot  consists  of  the  location  of  the 
critical  trajectories  in  the  y-z  plane  at  a  constant  value 
of  X  (distance  in  front  of  the  sampling  inlet)  equal  to  XDIN 
in  equation  (24).  In  this  case,  x  is  constant  because  the 
sampling  angle  is  0®  and  the  freestream  is  flowing  in  a 
direction  perpendicular  to  the  plane  at  x  =  XDIN. 

The  plot  in  figure  8  represents  the  area  in  the 
freestream,  Ao,  which  is  perpendicular  to  the  freestream 
resultant  velocity  for  the  case  when  the  sampling  angle  was 
30°.  In  this  case,  a  simple  plot  of  the  y-z  coordinates 
would  not  represent  the  actual  area  of  Ao  as  x  was  not 
constant  over  the  plane  of  interest.  To  solve  this  problem 
the  critical  trajectory  coordinates  were  converted  to  a 
local  coordinate  system  centered  around  the  coordinates  of 
the  initial  capture  point.  First  the  distance,  s,  between 
the  first  capture  point  and  the  critical  trajectory 
coordinate  on  each  of  the  13  rays  defined  by  a,  the 
incremental  angle,  was  calculated  using  the  distance 
formula,  equation  (22).  After  the  distances  between  the 
points  were  calculated,  the  x  and  y  coordinates  of  the 
critical  trajectories  determined  using  the  following 
relations : 

Y'  =  s  cosa  (26) 

Z  '  =  s  s  i  na 

The  local  coordinates  were  then  used  to  numerically 
calculate  the  area  using  the  FORTRAN  program,  AREA.F. 
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RESULTS 

Several  combinations  of  free  stream  velocity,  Uo,  inlet 
size,  di,  inlet  sampling  velocity,  Ui,  sampling  pitch  angle 
(misalignment  to  the  free  stream),  and  particle  aerodynamic 
diameter,  dae,  were  run  through  the  2D  and  3D  programs. 

Table  1  summarizes  the  input  variables  which  coincide  with 
input  parameters  used  by  Wiener  (1987). 


TABLE  1 

Model  Input  Parameters 


di 

Uo 

Ui 

dae 

theta 

cm _ _ 

cm/sec 

cm/sec 

micron 

_ degrees 

.32 

250 

125 

5.0 

0 

1.0 

500 

250 

40.0 

30 

1000 

500 

1000 

The  two  dimensional  program  was  run  for  96  combinations 
of  particle  size  (dae),  velocity  ratio  (Uo/Ui),  sampling 
angle  (6),  and  inlet  diameter  (di).  The  results  are 
summarized  in  the  appendix. 

A  linear  regression  of  the  BIEM  predicted  aspiration 
efficiency  on  the  empirical  aspiration  efficiency  calculated 
from  the  unified  aspiration  model,  equation  (2),  was 
performed  for  all  2D  data  collected,  for  all  0° 
combinations,  and  all  30°  combinations.  The  results  are 
summarized  in  table  2.  Figures  9,  10,  and  11  are  plots  of 


the  2D  aspiration  simulation  versus  the  unified  aspiration 
model,  equation  (2). 
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TABLE  2 

Two  Dimensional  Model  Linear 
Eegression  Results 


CATEGORY 

SLOPE 

Y  INT 

LCL  i  ' 

UCL  ^  ■' 

R2  (  3  ) 

All 

Data 

0.842 

0.074 

0.788 

0.895 

0.917 

All 

0° 

0 . 963 

0 . 077 

0 . 941 

0.986 

0 . 993 

All 

30° 

0.664 

0.118 

0.599 

0 . 729 

0.911 

(1) 

LCL  - 

95% 

lower  confidence 

limit  for 

s  1  ope 

(2) 

UCL  - 

95% 

upper  confidence 

limit  for 

s  1  ope 

(3)  R2  -  correlation  coefficient 

The  3D  results  are  presented  in  the  appendix  and 
summarized  in  figures  12,  13,  and  14.  Table  3  contains  the 
linear  regression  results  which  compare  the  numerical 
prediction  of  aspiration  efficiency  with  the  unified  model 
of  Hangal  and  Willeke,  1990. 


TABLE  3 

Three  Dimensional  Model  Linear 
Regression  Results 


CATEGORY _ 

SLOPE  _ 

Y  INT 

LCL  ^ 

UCL  ■: 

R  2  -  . 

All  Data 

1 . 062 

-0.062 

1 . 026 

1 . 099 

0 . 995 

All  0° 

1 . 059 

-0.026 

1 . 036 

1 .082 

0 . 999 

All  30° 

1 . 064 

-0 . 096 

0 . 980 

1  .  148 

0 . 991 

(1)  LCL  -  95%  lower  confidence  limit  for  slope 

(2)  UCL  -  95%  upper  confidence  limit  for  slope 

(3)  R2  -  correlation  coefficient 
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All  20  (0  and  30  degrees) 


Figure  9.  2D  results,  all  data.  Plot  of  the  predicted 
efficiency  developed  from  the  BIEM  and  Alenlus  theories 
against  the  empirical  unified  aspiration  efficiency  model  of 
Hangal  and  Hi 11  eke  (1990). 
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Isoaxiai  Case  (0  degrees) 


Figvire  10.  2D  results,  isoaxiai  data.  Plot  of  the 
predicted  efficiency  developed  from  the  BIEM  and  Alenius 
theories  against  the  empirical  unified  aspiration  efficiency 
model  of  Hangal  and  Willeke  (1990). 


0  2  4  e  8 

Unified  Model  Aspiration  Efficiency 


2  4  6  8 

Unified  Model  Aspiration  Efficiency 


Figure  12.  3D  results,  all  data.  Plot  of  the  predicted 
efficiency  developed  from  the  BIEM  and  Alenius  theories 
against  the  empirical  unified  aspiration  efficiency  model 
Hangal  and  Will eke  (1990). 
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3D  Non-lsoaxlal  Case  (30  degrees) 


Unified  Mode!  Aspiration  Efficiency 


Figure  14.  3D  results,  non-isoaxlal  data.  Plot  of  the 
predicted  efficiency  developed  from  the  BIEM  and  Alenlus 
theories  against  the  empirical  unified  aspiration  efficiency 
model  of  Hangal  and  Vllleke  (1990). 
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DISCUSSION 


The  two  dimensional  model  showed  good  agreement  with 
the  unified  aspiration  model  for  the  isoaxial  conditions. 
From  table  2,  the  slope  of  the  regression  line  for  the  0° 
sampling  angle  is  C.96.  A  slope  of  1.0  would  indicate 
perfect  agreement  with  the  unified  model .  The  largest  error 
in  the  2d  model  occurred  when  the  velocity  ratio  was  less 
than  1.  In  this  case,  the  inlet  velocity  was  2  to  4  times 
greater  than  the  wind  velocity  leading  to  large  changes  in 
air  velocity  and  direction  as  the  particles  approached  the 
inlet  face. 

The  nonisoaxial  2D  model  was  tested  for  a  30°  sampling 
angle  and  showed  poor  agreement  with  the  unified  model.  The 
percent  difference  between  the  unified  model  aspiration 
efficiency  and  the  predicted  aspiration  efficiency  was  as 
large  as  -71%.  The  slope  of  the  regression  line  for  this 
set  of  conditions  was  .66. 

Overall,  the  2D  model  fairly  represented  the  aspiration 
efficiency  of  aerosols  into  a  thin  edged  inlet.  The 
regression  equation,  slope  of  0.84,  was  biased  due  to  the 
poor  agreement  of  the  30°  model  with  the  unified  model. 

As  shown  in  table  3,  the  three  dimensional  solution 
performed  significantly  better  in  all  categories  when 
compared  to  the  2D  model.  The  slopes  of  the  regression 
lines  were  close  to  1.0  for  all  cases  indicating  that  the 
computer  model  accurately  calculates  aspiration  efficiency 
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values  when  compared  to  the  unified  aspiration  efficiency 
model.  Table  4  shows  a  comparison  of  the  2D  percent  error 
with  that  of  the  3D  model  for  the  same  combinations  of 
velocity  ratios,  sampling  angle  and  particle  diameter.  As 
can  be  seen,  in  most  cases,  the  3D  model  led  to  an 
improvement  in  the  aspiration  efficiency  value. 


TABLE  4 


Comparison  of 

2D 

Model  Error 

with  the 

3D 

Results 

uo 

UI 

THETA 

DP 

2D  ERR 

3D  e; 

1000 

125 

0 

5 

22% 

2% 

1000 

125 

0 

40 

4% 

5% 

250 

1000 

0 

5 

12% 

-9% 

250 

1000 

0 

40 

26% 

9% 

1000 

125 

30 

40 

-39% 

5% 

250 

1000 

30 

5 

-11% 

-11% 

250 

1000 

30 

40 

-3% 

-10% 

500 

500 

0 

5 

1% 

3% 

500 

500 

30 

5 

-10% 

0% 

500 

250 

0 

5 

11% 

6% 

500 

250 

30 

5 

-6% 

7% 

500 

125 

0 

5 

26% 

4% 

500 

125 

30 

5 

17% 

10% 

Some  error  was  introduced  into  the  above  comparison 
because  the  unified  model  of  aspiration  efficiency  was 
developed  for  circular,  thin  edged  inlets.  This  may  have 
some  effect  on  the  comparison  of  calculated  efficiency  with 
the  aspiration  efficiency  derived  from  the  computer.  This 
does  not,  however,  affect  the  model  calculation  of 
efficiency  because  the  actual  areas  using  the  square  inlet 
were  used  in  equation  (18). 
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Boundary  discretization,  particle  start  point, 
convergence  criteria,  and  time  step  (f)  were  the  most 
important  parameters  in  determining  accurate  solutions. 
Boundary  discretization  had  to  be  dense  enough  so  that 
calculations  of  internal  velocities  were  performed  at 
distances  greater  than  the  distance  between  nodes  on  the 
nearest  boundary  (Flynn  and  Miller,  1989).  This  was 
accomplished  by  developing  a  2D  boundary  consisting  of  75 
nodes  with  dense  discretization  near  the  inlet  face  and 
along  the  inlet  sides.  The  appendix  contains  a  sample  input 
file  for  the  2D  boundary  discretization.  With  this 
discretization,  few  aberrations  in  internal  velocity 
calculations  wers  seen.  This  was  not  the  case  for  the  3D 
solution. 

The  three  dimensional  model  boundary  discretization 
was  proportional  to  the  2D  discretization  with  the  sides  of 
the  triangular  elements  equal  in  length  to  the  distance 
between  nodes  in  the  2D  discretization.  This  led  to  an 
input  file  containing  367  nodes  defining  693  elements.  As 
described  by  Flynn  and  Miller,  the  triangular  elements  are 
defined  by  numbering  the  element  with  the  node  numbers 
corresponding  to  the  vertices  of  the  triangular  element  in  a 
clockwise  fashion  looking  out  of  the  domain.  This  was  done 
for  approximately  1/2  of  the  elements  starting  at  the  inlet 
face  and  along  the  walls  of  the  inlet  using  a  triangular 
element  generating  program  written  for  this  problem.  The 
remainder  of  the  element  input  file  was  fashioned  using 
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LOTUS  spreadsheets  because  of  the  complex  interactions  of 
common  nodes  along  the  edges  of  each  facet  of  the  cube  and 
at  its  corners. 

Even  with  this  dense  3D  boundary  discretization, 
internal  calculation  of  air  velocities  at  points  less  than 
0.05  cm  from  the  face  and  walls  of  the  inlet  were  extremely 
inaccurate  with  the  values  being  orders  of  magnitude  larger 
than  the  freestream  velocity.  This  led  to  large  errors  in 
the  prediction  of  particle  location  which  is  dependent  in 
part  on  accurate  calculation  of  air  velocity  at  the 
particle's  current  location.  To  avoid  this  erio^ ,  capture 
criteria  for  the  particle  was  extended  by  a  factor  of  0.05 
cm  in  front  of  the  inlet  face  and  around  its  edges.  This  in 
effect  enlarged  the  face  area,  Ai ,  to  1.11  cm^  for  a  square 
inlet  of  side  equal  to  1.0  cm.  This  enlarged  area  was  the 
area  used  in  the  aspiration  efficiency  equation  (18). 

The  effect  of  enlarging  this  area  was  to  alter  the 
actual  average  air  velocity  moving  into  the  "inlet".  Over 
the  face  of  the  inlet  the  air  velocity  was  specified  as  Ui . 
Over  the  thin  layer  defining  the  addition  to  the  actual 
inlet,  the  average  velocity  had  a  value  which  was  a  function 
of  the  free  stream  velocity  and  the  sampling  velocity. 
Because  of  this  there  is  some  error  in  the  specification  of 
Ui  in  equation  (164)  and  hence  Ea . 

Another  problem  with  the  3D  boundary  discretization 
involves  the  interface  of  a  surface  of  the  boundary  where 
the  velocity  is  specified  with  another  surface  where  the 
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velocity  is  also  specified,  but  is  different  than  that  of 
the  first  surface.  In  the  final  model,  this  occurred  only 
at  the  nodes  surrounding  the  edges  of  the  face  of  the 
sampling  inlet.  Here,  the  velocity  over  the  face  was 
specified  as  the  sampling  velocity,  Ui,  and  the  velocity  on 
the  edge  of  the  tube  and  along  its  length  was  specified  as  0 
cm/sec.  This  caused  a  discontinuity  of  velocity  at  the 
edges  of  the  sampling  inlet.  The  effect  of  the 
discontinuity  was  to  cause  inaccuracies  in  interpolation  of 
velocities  over  the  elements  involved.  This  interpolation 
error  was  minimized  by  inserting  a  very  thin  row  of  elements 
just  inside  the  edges  of  the  inlet.  For  the  final  model  the 
width  of  this  small  layer  was  0.001  cm. 

Particle  start  points  were  important  because  the 
aspiration  efficiency  calculation  requires  the  critical 
trajectory  start  points  to  be  in  the  free  stream.  It  was 
found  that  at  distances  of  6  inlet  diameters  from  the  inlet, 
local  air  velocity  and  direction  was  within  5%  of  the 
specified  free  stream  velocity.  This  agrees  with  the 
findings  of  Belyaev  and  Levin  (1972)  and  held  true  for  both 
2D  and  3D  models. 

When  determining  the  location  of  critical  trajectories, 
the  distance  between  the  last  capture  trajectory  start  point 
and  last  miss  trajectory  start  point  was  compared  to  the 
convergence  criteria.  Once  the  distance  reached  a  preset 
criteria,  the  points  were  averaged  to  give  the  start  point 
of  the  critical  trajectory.  The  convergence  criteria  for 
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the  2D  model  was  set  at  0.01  cm  after  iterations  at  criteria 
smaller  than  this  caused  a  less  than  1%  change  in  the 
calculated  efficiency.  The  convergence  criteria  for  the  3D 
model  led  to  good  agreement  with  the  unified  model  when  set 
at  0.05  cm.. 

The  time  step  was  the  most  critical  factor  in 
developing  accurate  particle  trajectories.  Time  steps  that 
were  too  large  resulted  in  predictions  of  trajectories  that 
were  too  coarse  to  accurately  predict  the  particle  movement. 
Unnecessarily  small  time  steps  caused  excessive  computing 
time.  Through  a  series  of  program  tests,  0.0003  seconds  was 
chosen  as  the  time  step  for  the  2D  model  with  a  further 
reduction  by  a  factor  of  10  once  the  particle  came  within  1 
inlet  diameter  of  the  sampling  inlet  to  account  for  the 
large  velocity  changes  near  the  face  of  the  inlet. 

The  time  step  and  the  convergence  criteria  for  the  3D 
model  were  not  as  sensitive  in  the  determination  of  the 
aspiration  efficiency  as  they  were  in  the  2D  case.  Since 
the  3D  model  more  realistically  represented  the  physics  of 
the  problem,  accurate  solutions  were  obtained  for  time  steps 
ranging  from  0.001  to  0.0003  seconds  with  the  time  step 
being  halved  when  the  particle  was  within  1  inlet  diameter 
of  the  inlet  face. 

Because  the  size  of  the  time  step  indirectly  predicts 
the  distance  the  particle  will  move;  positions  predicted 
when  the  air  is  turning  and  accelerating  will  not  lie  on  the 
the  actual  trajectory.  This  is  clearly  seen  in  the  3D  model 
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for  velocity  ratios  of  2  and  4  where  the  estimate  of  Ea  was 
consistently  larger  than  the  unified  model  prediction.  In 
this  case,  the  air  is  approaching  the  inlet  faster  than  the 
inlet  is  sampling  and  must  diverge  around  the  inlet. 

Particle  trajectories  will  also  diverge  but  the  predicted 
position  will  lie  slightly  inside  of  the  actual  trajectory. 
This  causes  an  overestimate  of  Ao  leading  to  an  increase  in 
the  prediction  of  Ea. 

When  running  the  model  under  certain  conditions,  the 
time  step  and  capture  criteria  became  very  critical  in 
determining  whether  a  trajectory  resulted  in  a  miss  or 
capture.  In  the  three  dimencional  model,  when  Uo  =  1000,  Ui 
=  125  cm/sec,  and  0  =  30°,  the  best  estimate  of  aspiration 
efficiency  was  off  by  39%  when  compared  to  the  unified 
model.  Under  these  conditions  the  air  flow  is  approaching 
the  inlet  at  8  times  the  sampling  rate  causing  significant 
divergence  of  the  air  streamlines  around  the  face  of  the 
inlet.  This  situation  combined  with  the  sampling  angle 
results  in  particle  trajectories  which  move  virtually 
straight  up  in  the  positive  z  direction  at  locations  within 
0.1  cm  of  the  inlet  face  when  the  particle  is  approaching 
the  inlet  from  above  the  centerline.  The  air  velocity  in 
this  area  is  also  moving  very  fast.  The  combination  of  high 
velocity  particle  and  steep  angle  of  approach  toward  the 
sample  inlet  can  easily  lead  to  a  predicted  particle 
location  which  jumps  across  the  numerical  edges  of  the 
sampling  inlet.  A  miss  trajectory  is  recorded  ^or  this 
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situation  when  in  actuality  a  capture  should  have  occurred. 
This  causes  an  underestimation  of  Ao  leading  to  a  prediction 
of  Ea  which  is  too  small.  This  situation  occurred  for  the 
last  four  rays  around  the  start  point  semicircle  and  is  the 
reason  for  the  poor  performance  of  the  model  under  these 
conditions.  This  situation  can  be  eliminated  by  reducing 
the  size  of  the  time  step  but  would  lead  to  a  large  increase 
in  computing.  The  same  problem  was  seen  in  the  2D  model  and 
prevented  identification  of  an  upper  critical  trajectory. 

By  specifiying  an  incremental  angle  of  15°  in  the  3D 
model,  the  program  was  limited  to  locating  13  critical 
trajectories.  The  effect  of  decreasing  a  has  not  been 
analyzed  but  would  lead  to  a  finer  outline  of  the  critical 
trajectory  stream  tube  and  a  more  accurate  calculation  of 
Ao .  One  anomaly  in  the  output  was  seen  when  the  ray 
specified  conincided  with  the  corners  of  the  inlet. 

Critical  trajectories  associated  with  these  rays  tended  to 
lie  outside  range  of  critical  trajectories  on  either  side  of 
this  point.  In  one  program  test,  eliminating  these  rays 
from  the  program  improved  the  model  result  from  a  26%  to  a 
4%  error. 

The  results  presented  in  this  paper  show  that  the 
solutions  agree  with  the  unified  empirical  model  of  Hangal 
and  Willeke  (1990).  The  advantage  of  this  model  is  that  it 
can  be  expanded  to  situations  beyond  that  investigated. 

These  situations  include  modeling  a  yaw  angle,  multiple 
inlets,  and  inlets  of  different  geometry.  Modeling  a  yaw 
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angle  as  opposed  to  a  pitch  angle  would  require  only  a  shift 
in  the  direction  of  the  gravity  component  in  the  currently 
available  boundary  file.  Modeling  of  multiple  inlets  and 
differing  inlet  geometries  could  be  accomplished  by 
reworking  boundary  files  to  represent  the  new  inlet 
configurations . 
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CONCLUSION 

The  boundary  integral  equation  method  can  be  used  to 
accurately  describe  the  air  flow  field  into  a  thin  edged 
inlet.  Once  the  flow  is  described,  Alenius's  (1989) 
analytical  solution  to  the  particle  equations  of  motion  can 
be  used  to  locate  critical  aerosol  trajectories  and 
accurately  predict  aspiration  efficiency.  A  simple  two 
dimensional  model  described  here  can  accurately  predict 
aspiration  efficiency  for  small  particles  when  the  velocity 
ratio  is  close  to  1 .  In  this  case,  the  critical 
trajectories  are  symmetric  about  the  center  of  the  inlet. 
When  an  angle  of  misalignment  is  introduced  and/or  the 
velocity  ratios  differ  significantly  from  1.0,  a  three 
dimensional  model  is  required  to  describe  critical 
trajectories  and  thus  aspiration  efficiency. 

The  3D  model  can  accurately  predict  aspiration 
efficiencies  which  are  in  good  agreement  with  the  unified 
aspiration  efficiency  model  for  all  forward  sampling  angles 
proposed  by  Hangal  and  Willeke  (1990).  Further  refinements 
to  the  model  should  improve  the  accuracy  of  aspiration 
efficiency  estimates. 

Finally,  this  model  has  the  capability  of  predicting 
aspiration  efficiencies  for  inlets  of  differing  shapes  and 
orientations  including  multiple  inlets,  yaw  angles,  and 
inlets  of  differing  geometry.  By  reworking  the  boundary 
input  files,  inlets  of  any  geometry,  location,  sampling 
angle  and  size  can  be  studied  using  these  programs. 
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FUTURE  WORK 

Successful  calculation  of  particle  trajectories  depends 
on  several  adjustable  variables  in  the  program  including 
boundary  discretization,  calculation  time  step,  convergence 
criteria,  and  capture  criteria.  While  this  paper 
demonstrates  that  the  method  can  accurately  track  particles 
for  given  conditions,  the  above  variable  have  not  been 
tested  to  determine  their  optimum  values.  A  study  of  these 
variables  will  help  determine  the  contribution  of  each  to 
the  overall  error  in  the  program. 

Boundary  discretization  of  the  3D  domain  required  to 
solve  the  BIEM  problem  can  be  difficult  due  to  the  number  of 
triangular  elements  required.  Also,  when  constructing  a 
closed  domain,  a  number  of  common  nodes  appear  along  the 
edges  of  each  face  of  the  boundary  and  at  the  corners.  This 
can  lead  to  cases  where  one  node  is  contained  in  up  to  6 
distinct  triangular  elements.  Boundary  files  for  this 
project  were  manually  prepared  for  about  50%  of  the  domain. 
Full  automation  of  this  process  would  greatly  increase  the 
applicability  of  this  method. 

Finer  discretization  at  the  face  of  the  inlet  and  along 
the  walls  of  the  inlet  will  decrease  the  distance  a  particle 
can  come  within  the  inlet  and  inlet  sides  before  internal 
air  velocity  calculations  no  longer  represent  the 
conditions.  In  this  project,  that  distance  was  0.05  cm. 

This  will  lead  to  a  significant  increase  in  the  number  of 
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nodes  and  trian^jular  elements  required  to  describe  the 
boundary  and  may  increase  the  size  of  the  boundary  file  by 
as  much  as  50%  and  subsequently  increase  computing  time. 
Finer  discretization  may  be  avoided  by  adding  an 
interpolation  scheme  to  the  program.  This  subroutine  would 
be  invoked  when  the  particle  entered  the  volume  around  the 
face  of  the  inlet  where  internal  velocity  calculations 
showed  large  error.  While  this  will  cause  some  error  in  the 
prediction  of  internal  velocity,  it  is  expected  it  will  be 
less  than  that  of  assuming  a  captured  trajectory  has 
occurred  when  the  particle  enters  this  area. 

Rediscretization  of  the  inlet  and  walls  to  more  closely 
represent  a  circular  inlet  may  also  lead  to  a  better 
comparison  of  the  model  with  the  unified  aspiration 
efficiency  equation.  This  could  be  done  by  modelling  the 
inlet  as  an  equal  area  hexagon  or  octagon.  Automation  of 
the  boundary  file  would  be  paramount  in  this  case  as  the 
discretization  would  become  extremely  tedious. 

It  may  be  necessary  to  add  a  routine  to  the  program  to 
automatically  decrease  the  time  step  if  the  distance  between 
successive  points  becomes  large  or  the  particle  jumps  across 
the  boundaries  of  the  inlet.  Alternatively,  the  path  of  the 
particle  can  be  analyzed  to  see  if  it  crossed  the  inlet  and 
a  routine  added  to  the  capture  criteria  which  records  a 
capture  when  this  situation  occurs. 

The  boundary  input  files  and  start  point  routines  were 
designed  for  positive  sampling  angles  only.  The  program  is 
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not  general  enough  to  support  negative  (downward  pointing) 
inlet  orientations.  However,  by  simply  changing  the  sign  on 
the  gravity  vector  and  representing  the  negative  angle  as  a 
positive  angle,  predictions  of  aspiration  efficiency  may  be 
obtained . 

Combination  of  the  3D  boundary  solution  program  with 
the  particle  tracking  routine  and  the  area  integration 
routine  would  further  streamline  the  efficiency 
calculations.  Currently,  3  FORTRAN  programs  are  used  to 
calculate  aspiration  efficiency. 
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APPENDIX 
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Calculation  of  Particle  Trajectories : 


Particle  transport  is  influenced  by  a  gravity  force  and 
a  drag  force  due  to  the  friction  of  the  svirroundlng  air.  By 
Newton's  second  law,  particle  mass  times  acceleration  is 
equal  to  the  sum  of  the  forces  on  the  particle, 

€F  =  ma  (1) 


where  F  is  the  force,  a  is  the  acceleration  and  bold  face 
quantities  represent  vectors. 


The  frictional  force  can  be  expressed  by  the  following 
equation: 


Fd  =  -  Cd(Re)  Ap  pi  vr  vr  (2) 

CC  2 


where 


Cd(Re)  is  the  dimensionless  friction  factor 
dependent  on  the  Reynolds  number 

CC  is  the  Cunningham  slip  correction 

Ap  is  the  projected  area  of  the  particle 

pi  is  the  density  of  the  medium  (air) 

vr  is  the  magnitude  of  the  particle  relative 
velocity 
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vr  is  the  particle  velocity  relative  to  the 
air. 


Cd(Re)  is  dependent  on  the  particle  Reynold's  number. 

Re  =  vr  dp  pi  (3) 

Ml 

where 


dp  is  the  particle  diameter  (assumed 
spherical ) 

pi  is  the  viscosity  of  the  air. 


Several  expressions  for  Cd(Re)  can  be  found  in  the 
literature  (Relst,  1964)  with  most  yielding  no  major 
deviation  (Alenlus,  1967).  For  this  project  the  expressions 
used  by  Alenius  are  Incorporated  into  a  simplifying 
expression: 


Cd(Re)=24/Re  Q(Re) 


where  Q(Re)  is  defined  as: 


(4) 


Q(Re)  =  1  for  0  <  Re  <=  0.5 
Q(Re)  =  1  +  0.15  RoO  t-e?  for  0.5  <  Re  <=  800 
Q(Re)  =  0.44  Re/24  for  600  <  Re  <  20000. 


By  combining  equations  (3)  and  (4)  into  equation  (2) 
and  Introducing  the  expression  for  particle  relaxation  time, 
T,  a  simplified  expression  for  the  drag  force  becomes: 
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where 


Fd/m  =  Q(Re)/T  vr 


(5) 


T  =  dp»  pi  CC  (6) 

18  Ml 

From  equations  (1)  and  (5)  the  particle  motion  can  be 
described  as: 

d^r/dti  -  [Q(Re)/T(u-(dr/dt))]-g  =  0  (7) 

where 

r  is  the  particle  position  vector 

g  is  the  acceleration  due  to  gravity 

u  is  the  air  velocity  at  the  point  where 
the  particle  is  located. 

Equation  (7)  is  analytically  solvable  if  the  factor 
Q(Re)  is  constant  and  the  air  velocity  u  at  the  points  where 
the  particle  is  successively  located  can  be  described  as  a 
linear  function  of  the  time  that  has  elapsed  until  the 
particle  reached  the  point.  This  condition  can  be 
approximately  met  it  the  motion  is  calculated  in 
sufficiently  small  time  steps. 

If  we  introduce  an  expression  for  the  linear  behavior 
of  the  air  velocity  over  the  time  interval  being  studied. 


then: 
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u(t)  =  (u(t)-u(t-t' )/t'+b  (8) 

where 

t  is  time 

t'  is  the  length  of  the  time  step 
b  is  some  initial  value  of  velocity 

For  ease  of  writing  we  define  the  first  part  of  the  above 
expression  as  ua: 

ua  =  (u(t)-u(t-t' )/t ’ )  (9) 

Now  if  we  let  A=Q(Re)/T  equation  (7)  can  be  rewritten 

as : 

dv/dt  -  A(uat+b)  +  Av  -  g  =  0 
recognizing  that  d*r/dt*  =  dv/dt, 

or,  dv/dt  +  Av  =  A(uat+b)-  g  (10) 

If  equation  (10)  is  multiplied  through  by  e*t  the  left 
hand  side  of  the  equation  will  be  a  derivative. 

eAt (dv/dt)  +  e*t(Av)  =(e» t ) A(uat+b)-  g(eAt) 

Integrating  both  sides  of  the  above  equation  from  tl  to  t2 
gives : 

V( t 2 ) e* t 2-V{ t 1 ) e* t 1  ■  [(uat2)e*t2  -  (ua/A)e*t2 

+(b)e*t2  +  (g/A)eAt2] 

-  [(uat2)e*ti  -  (ua/A)eAti 


+(b)e*ti  +  (g/A)eAti] 


(11) 
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Adding  V(ti)eAti  to  both  side  of  the  equation,  dividing 
through  by  e*t2  and  simplifying  gives  an  expression  for 

V(  1 2 )  : 

V(t2)  =  uat2  -  ua/A  +  b  +  g/A 

-  {(uatl)  “  ua/A  +  b  +  g/A  ~  (12) 

If  we  define  a  simplifying  term, 

Ub  =  u(t)-(Ua-g)T/Q(Re) 

and  define  tl  as  t,  and  t2  as  the  time  at  the  end  of  the 
times  step,  t'  (t2  =  t  +  t*),  substitute  the  air  velocity  at 
time  t  for  b,  and  replace  A  with  Q(Re)/T,  then  the  solution 
is: 

v(t+t')  =  u«t '+Ub-(Ub-v(t)  )exp- [Q(Re)t ' /'^]  (13) 

Integrating  the  equation  (13)  again  will  give  the  particle 
position: 

r(t+t*)  =  r(t)+u./2t ' 2+Ubt’ 

-(Ub-v(t))T/Q(Re) (l-expQ(Re)t'/T)  (14) 

Equation  (13)  and  (14)  are  analytical  solutions  to  the 
particle  equations  of  motion  and  are  used  in  both  the  2D  and 
3D  particle  tracking  subroutines. 
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2D  Boundary  Discretization  Sample: 


Potential  specified  on  top  and  rear  sides. 
Velocity  specified  over  front,  bottom,  inlet  face 
and  sides. 


NODE 

DIN 

THETA 

UI 

UO 

UOX 

UOY 

X 

1 

30 

1000 

-250 

-216.5064 

125 

Y 

cm 

degrees 

cm/sec 

cm/sec 

PHI 

(inlet  diameter) 
(nozzle  angle) 
(sampling  velocity) 
(wind  velocity) 

(x  component  of  UO) 

(y  component) 

DPHI  ID  ICT 

1 

20 

15 

0 

-216.5064 

3 

0 

2 

20 

10 

0 

-216.5064 

3 

0 

3 

20 

5 

0 

-216.5064 

3 

0 

4 

20 

0 

0 

-216.5064 

3 

0 

5 

20 

-5 

0 

-216.5064 

3 

0 

6 

20 

-10 

0 

-216.5064 

3 

0 

7 

20 

-15 

0 

-216.5064 

3 

0 

8 

20 

-20 

0 

-216.5064 

3 

0 

9 

19.99 

-20 

0 

-125 

3 

0 

10 

15 

-20 

0 

-125 

3 

0 

11 

10 

-20 

0 

-125 

3 

0 

12 

5 

-20 

0 

-125 

3 

0 

13 

0 

-20 

0 

-125 

3 

0 

14 

-5 

-20 

0 

-125 

3 

0 

15 

-10 

-20 

0 

-125 

3 

0 

16 

-15 

-20 

0 

-125 

3 

0 

17 

-20 

-20 

1830. 

13 

0 

2 

1 

18 

-20 

-15 

2455. 

13 

0 

2 

0 

19 

-20 

-10 

3080. 

13 

0 

2 

0 

20 

-20 

-5 

3705. 

13 

0 

2 

0 

21 

-20 

-0.5 

4267. 

63 

0 

2 

2 

22 

-17.5 

-0.5 

0 

0 

3 

0 

23 

-15 

-0.5 

0 

0 

3 

0 

24 

-12.5 

-0.5 

0 

0 

3 

0 

25 

-10 

-0.5 

0 

0 

3 

0 

26 

-7.5 

-0.5 

0 

0 

3 

0 

27 

-5 

-0.5 

0 

0 

3 

0 

28 

-4.5 

-0.5 

0 

0 

3 

0 

29 

-4 

-0.5 

0 

0 

3 

0 

30 

-3.5 

-0.5 

0 

0 

3 

0 

31 

-3 

-0.5 

0 

0 

3 

0 

32 

-2.5 

-0.5 

0 

0 

3 

0 

33 

-2 

-0.5 

0 

0 

3 

0 

34 

-1.5 

-0.5 

0 

0 

3 

0 

35 

-1 

-0.5 

0 

0 

3 

0 

60 


36 

-0.5 

-0.5 

0 

37 

0 

-0.5 

0 

3B 

0 

-0.499 

0 

39 

0 

-0.375 

0 

40 

0 

-0.25 

0 

41 

0 

-0.125 

0 

42 

0 

0 

0 

43 

0 

0.125 

0 

44 

0 

0.25 

0 

45 

0 

0.375 

0 

46 

0 

0.499 

0 

47 

0 

0.5 

0 

48 

-0.5 

0.5 

0 

49 

-1 

0.5 

0 

50 

-1.5 

0.5 

0 

51 

-2 

0.5 

0 

52 

-2.5 

0.5 

0 

53 

-3 

0.5 

0 

54 

-3.5 

0.5 

0 

55 

-4 

0.5 

0 

56 

-4.5 

0.5 

0 

57 

-5 

0.5 

0 

58 

-7.5 

0.5 

0 

59 

-10 

0.5 

0 

60 

-12.5 

0.5 

0 

61 

-15 

0.5 

0 

62 

-17.5 

0.5 

0 

63 

-20 

0.5 

4392.63 

64 

-20 

5 

4955.13 

65 

-20 

10 

5580.13 

66 

-20 

15 

6205.13 

67 

-20 

20 

6830.13 

68 

-15 

20 

5747.6 

69 

-10 

20 

4665.06 

70 

-5 

20 

3582.53 

71 

0 

20 

2500 

72 

5 

20 

1417.47 

73 

10 

20 

334.94 

74 

15 

20 

-747 . 6 

75 

20 

20 

0 

0 

0 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

1000 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-216.5064 


3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

3 

2 

2 

2 

2 

2 

2 

2 

2 

2 


0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 


d  r»  CO 
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3D  Boundary  Discretization  Sample; 

INPUT  FILE  FOR  3D  CASE,  367  NODE,  693  ELEMENTS 

Potential  specified  at  all  faces  except  at  boudary  of  tube 
and  at  the  face.  f>ur  corners  specified  with  potential. 


DIN 

1  cm 

(inlet  diameter) 

THETA 

30  degrees 

(nozzle  angle) 

UI 

1000  cm/sec 

(sampling  velocity) 

UO 

-250  cm/sec 

(wind  velocity) 

UOX 

-216.5064 

(X  component  of  UO) 

UOY 

0 

(y  component) 

UOZ 

125 

(z  component) 

62 


NODE  X  Y 


1 

0 

-0.5 

2 

0 

-0.5 

3 

0 

-0.5 

4 

0 

-0.5 

5 

0 

-0.5 

6 

0 

-0.5 

7 

0 

-0.5 

8 

0 

-0.5 

9 

0 

-0.5 

10 

0 

-0.5 

1 1 

0 

-0.5 

1  2 

0 

-0.499 

13 

0 

-0.499 

14 

0 

-0.499 

15 

0 

-0.499 

16 

0 

-0.499 

17 

0 

-0.499 

18 

0 

-0.499 

19 

0 

-0.499 

20 

0 

-0.499 

21 

0 

-0.499 

22 

0 

-0.499 

23 

0 

-0.375 

24 

0 

-0.375 

25 

0 

-0.375 

26 

0 

-0.375 

27 

0 

-0.375 

28 

0 

-0.375 

29 

0 

-0.375 

30 

0 

-0.375 

31 

0 

-0.375 

32 

0 

-0.375 

33 

0 

-0.375 

34 

0 

-0.25 

35 

0 

-0.25 

36 

0 

-0.25 

37 

0 

-0.25 

38 

0 

-0 . 25 

39 

0 

-0.25 

40 

0 

-0.25 

41 

0 

-0.25 

42 

0 

-0 . 25 

43 

0 

-0.25 

44 

0 

-0.25 

45 

0 

-0.125 

46 

0 

-0.125 

47 

0 

-0.125 

48 

0 

-0.125 

49 

0 

-0.125 

50 

0 

-0.125 

51 

0 

-0.125 

52 

0 

-0.125 

53 

0 

-0.125 

2 

PHI 

DPHI 

ID 

ICT 

-0 

0 

0 

2 

0 

-0.499 

0 

0 

2 

0 

-0.375 

0 

0 

2 

0 

-0.25 

0 

0 

2 

0 

-0.125 

0 

0 

2 

0 

0 

0 

0 

2 

0 

0.125 

0 

0 

2 

0 

0.25 

0 

0 

2 

0 

0.375 

0 

0 

2 

0 

0.499 

0 

0 

2 

0 

0 . 5 

0 

0 

2 

0 

-0.5 

0 

0 

2 

0 

-0.499 

0 

1000 

2 

0 

-0.375 

0 

1000 

2 

0 

-0.25 

0 

1000 

2 

0 

-0.125 

0 

1000 

2 

0 

0 

0 

1000 

2 

0 

0.125 

0 

1000 

2 

0 

0.25 

0 

1000 

2 

0 

0.375 

0 

1000 

2 

0 

0.499 

0 

1000 

2 

0 

0.5 

0 

0 

2 

0 

-0.5 

0 

0 

2 

0 

-0.499 

0 

1000 

2 

0 

-0.375 

0 

1000 

2 

0 

-0.25 

0 

1000 

2 

0 

-0.125 

0 

1000 

2 

0 

0 

0 

1000 

2 

0 

0. 125 

0 

1000 

2 

0 

0.25 

0 

1000 

2 

0 

0.375 

0 

1000 

2 

0 

0.499 

0 

1000 

2 

0 

0.5 

0 

0 

2 

0 

-0.5 

0 

0 

2 

0 

-0.499 

0 

1000 

2 

0 

-0.375 

0 

1000 

2 

0 

-0.25 

0 

1000 

2 

0 

-0.125 

0 

1000 

2 

0 

0 

0 

1000 

2 

0 

0.125 

0 

1000 

2 

0 

0.25 

0 

1000 

2 

0 

0.375 

0 

1000 

2 

0 

0.499 

0 

1000 

2 

0 

0.5 

0 

0 

2 

0 

-0 . 5 

0 

0 

2 

0 

-0.499 

0 

1000 

2 

0 

-0.375 

0 

1000 

2 

0 

-0.25 

0 

1000 

2 

0 

-0.125 

0 

1000 

2 

0 

0 

0 

1000 

2 

0 

0.125 

0 

1000 

2 

0 

0.25 

0 

1000 

2 

0 

0.375 

0 

1000 

2 

0 

63 


54 

0 

-0.125 

0.499 

0 

1000 

2 

0 

55 

0 

-0.125 

0.5 

0 

0 

2 

0 

56 

0 

0 

-0.5 

0 

0 

2 

0 

57 

0 

0 

-0.499 

0 

1000 

2 

0 

58 

0 

0 

-0.375 

0 

1000 

2 

0 

59 

0 

0 

-0.25 

0 

1000 

2 

0 

60 

0 

0 

-0.125 

0 

1000 

2 

0 

61 

0 

0 

0 

0 

1000 

2 

0 

62 

0 

0 

0.125 

0 

1000 

2 

0 

63 

0 

0 

0.25 

0 

1000 

2 

0 

64 

0 

0 

0.375 

0 

1000 

2 

0 

65 

0 

0 

0.499 

0 

1000 

2 

0 

66 

0 

0 

0.5 

0 

0 

2 

0 

67 

0 

0.125 

-0.5 

0 

0 

2 

0 

68 

0 

0.125 

-0.499 

0 

1000 

2 

0 

69 

0 

0.125 

-0.375 

0 

1000 

2 

0 

70 

0 

0.125 

-0.25 

0 

1000 

2 

0 

71 

0 

0-125 

-0.125 

0 

1000 

2 

0 

72 

0 

0.125 

0 

0 

1000 

2 

0 

73 

0 

0.125 

0.125 

0 

1000 

2 

0 

74 

0 

0.125 

0.25 

0 

1000 

2 

0 

75 

0 

0.125 

0.375 

0 

1000 

2 

0 

76 

0 

0.125 

0.499 

0 

1000 

2 

0 

77 

0 

0.125 

0.5 

0 

0 

2 

0 

78 

0 

0 . 25 

-0.5 

0 

0 

2 

0 

79 

0 

0.25 

-0.499 

0 

1000 

2 

0 

80 

0 

0.25 

-0.375 

0 

1000 

2 

0 

81 

0 

0.25 

-0.25 

0 

1000 

2 

0 

82 

0 

0.25 

-0.125 

0 

1000 

2 

0 

83 

0 

0.25 

0 

0 

1000 

2 

0 

84 

0 

0.25 

0.125 

0 

1000 

2 

0 

85 

0 

0.25 

0.25 

0 

1000 

2 

0 

86 

0 

0.25 

0.375 

0 

1000 

n 

0 

87 

0 

0.25 

0.499 

0 

1000 

2 

0 

88 

0 

0.25 

0.5 

0 

0 

2 

0 

89 

0 

0.375 

-0.5 

0 

0 

2 

0 

90 

0 

0.375 

-0.499 

0 

1000 

2 

0 

91 

0 

0.375 

-0.375 

c 

1000 

2 

0 

92 

0 

0.375 

-0 . 25 

0 

1000 

2 

0 

93 

0 

0.375 

-0.125 

0 

1000 

2 

0 

94 

0 

0.375 

0 

0 

1000 

2 

0 

95 

0 

0.375 

0.125 

0 

1000 

2 

0 

96 

0 

0.375 

0.25 

0 

1000 

2 

0 

97 

0 

0.375 

0.375 

0 

1000 

2 

0 

98 

0 

0.375 

0.499 

0 

1000 

2 

0 

99 

0 

0.375 

0.5 

0 

0 

2 

0 

100 

0 

0.499 

-0.5 

0 

0 

2 

0 

101 

0 

0.499 

-0.499 

0 

1000 

2 

0 

102 

0 

0 . 499 

-0.375 

0 

1000 

2 

0 

103 

0 

0.499 

-0.25 

0 

1000 

2 

0 

104 

0 

0.499 

-0.125 

0 

1000 

2 

0 

105 

0 

0.499 

0 

0 

1000 

2 

0 

106 

0 

•  0.499 

0.125 

0 

1000 

2 

0 

107 

0 

0.499 

0.25 

0 

1000 

2 

0 

64 


108 

0 

0.499 

0.375 

0 

1000 

2 

0 

109 

0 

0.499 

0.499 

0 

1000 

2 

0 

1  1  0 

0 

0.499 

0 . 5 

0 

0 

2 

0 

1  1  1 

0 

0 . 5 

-0.5 

0 

0 

2 

0 

1  1  2 

0 

0.5 

-0.499 

0 

0 

2 

0 

113 

0 

0.5 

-0.375 

0 

0 

2 

0 

114 

0 

0 . 5 

-0  .  25 

0 

0 

2 

0 

115 

0 

0.5 

-0-125 

0 

0 

2 

0 

116 

0 

0 . 5 

0 

0 

0 

2 

0 

117 

0 

0.5 

0.125 

0 

0 

2 

0 

118 

0 

0 . 5 

0.25 

0 

0 

2 

0 

119 

0 

0.5 

0.375 

0 

0 

2 

0 

120 

0 

0.5 

0.499 

0 

0 

2 

0 

121 

0 

0.5 

0.5 

0 

0 

2 

0 

122 

-0 . 5 

-0.5 

-0.5 

0 

0 

2 

0 

123 

-1 

-0.5 

-0.5 

0 

0 

2 

0 

124 

-1  .  5 

-0.5 

-0.5 

0 

0 

2 

0 

125 

-2 

-0.5 

-0 . 5 

0 

0 

2 

0 

126 

-2.5 

-0.5 

-0 . 5 

0 

0 

2 

0 

127 

-3 

-0 . 5 

-0.5 

0 

0 

2 

0 

128 

-3 . 5 

-0 . 5 

-0.5 

0 

0 

2 

0 

129 

-4 

-0.5 

-0.5 

0 

0 

2 

0 

130 

-4.5 

-0.5 

-0 . 5 

0 

0 

2 

0 

131 

-5 

-0 . 5 

-0.5 

0 

0 

2 

0 

132 

-0.5 

0 

-0 . 5 

0 

0 

2 

0 

133 

-1 

0 

-0.5 

0 

0 

2 

0 

134 

-1  .5 

0 

-0.5 

0 

0 

2 

0 

135 

-2 

0 

-0.5 

0 

0 

2 

0 

136 

-2 . 5 

0 

-0.5 

0 

0 

2 

0 

137 

-3 

0 

-0.5 

0 

0 

2 

0 

138 

-3 . 5 

0 

-0 . 5 

0 

0 

2 

0 

139 

-4 

0 

-0 . 5 

0 

0 

2 

0 

140 

-4.5 

0 

-0 . 5 

0 

0 

2 

0 

141 

-5 

0 

-0.5 

0 

0 

2 

0 

142 

-0.5 

0.5 

-0 . 5 

0 

0 

2 

0 

143 

-1 

0.5 

-0.5 

0 

0 

2 

0 

144 

-1  .  5 

0.5 

-0 . 5 

0 

0 

2 

0 

145 

-2 

0.5 

-0.5 

0 

0 

2 

0 

146 

-2 . 5 

0.5 

-0.5 

0 

0 

2 

0 

147 

-3 

0.5 

-0.5 

0 

0 

2 

0 

148 

-3 . 5 

0 . 5 

-0 . 5 

0 

0 

2 

0 

149 

-4 

0 . 5 

-0 . 5 

0 

0 

2 

0 

150 

-4.5 

0.5 

-0.5 

0 

0 

2 

0 

151 

-5 

0 . 5 

-0 . 5 

0 

0 

c 

0 

152 

-0 . 5 

0 . 5 

0 

0 

0 

2 

0 

153 

-1 

0.5 

0 

0 

0 

2 

0 

1  54 

-1.5 

0.5 

0 

0 

0 

2 

0 

155 

-2 

0 . 5 

0 

0 

0 

2 

0 

156 

-2 . 5 

0 . 5 

0 

0 

0 

2 

0 

157 

-3 

0.5 

0 

0 

0 

2 

0 

158 

-3 . 5 

0.5 

0 

0 

0 

2 

0 

159 

-4 

0.5 

0 

0 

0 

2 

0 

160 

-4 . 5 

0 . 5 

0 

0 

0 

2 

0 

161 

-5 

0.5 

0 

0 

0 

2 

0 

65 


162 

-0 . 5 

0.5 

0.5 

0 

0 

2 

0 

163 

-1 

0.5 

0.5 

0 

0 

2 

0 

164 

-1  .  5 

0.5 
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Two  Dimensional  Model  Output 


Time  Step  =  0.003  seconds 
Convergence  Criteria  =  0.001 


MODEL  UNIFIED 


UO.cm/sec 

Ui .cm/sec 

Theta  d . 

um 

D .  cm 

Eff . 

MODEL  EFF 

250 

125 

0 

5 

1  .  00 

1.17 

1  . 0420 

250 

125 

0 

40 

1.00 

1.71 

1  .7372 

500 

125 

0 

5 

1  .00 

1  .55 

1 . 2268 

500 

125 

0 

40 

1  .  00 

3.40 

3.5187 

1000 

125 

0 

5 

1  .  00 

2.39 

1  .  9534 

1000 

125 

0 

40 

1 . 00 

7 . 04 

7.3689 

250 

250 

0 

5 

1  .  00 

0.99 

1 . 0000 

250 

250 

0 

40 

1  .  00 

1  .  00 

1  .0000 

500 

250 

0 

5 

1  .  00 

1  .  20 

1 . 0806 

500 

250 

0 

40 

1  .00 

1.81 

1 . 8487 

1000 

250 

0 

5 

1  .  00 

1  .  69 

1.4217 

1000 

250 

0 

40 

1  .  00 

3.58 

3 . 7384 

250 

500 

0 

5 

1  .  00 

0 .89 

0.9710 

250 

500 

0 

40 

1  .  00 

0 . 64 

0.6013 

500 

500 

0 

5 

1  .  00 

0.99 

1 . 0000 

500 

500 

0 

40 

1  .  00 

1  .  00 

1 . 0000 

1000 

500 

0 

5 

1 . 00 

1  .  26 

1 . 1492 

1000 

500 

0 

40 

1  .  00 

1  .  87 

1.9182 

250 

1000 

0 

5 

1  .  00 

0.82 

0.9412 

250 

1000 

0 

40 

1  .00 

0.46 

0.3665 

500 

1000 

0 

5 

1.00 

0.86 

0 . 9453 

500 

1000 

0 

40 

1  .  00 

0.60 

0.5564 

1000 

1000 

0 

5 

1  .  00 

0 . 99 

1 . 0000 

1  000 

1000 

0 

40 

1  .  00 

1  .  00 

1 . 0000 

250 

125 

0 

5 

0.32 

1  .  25 

1.1205 

250 

125 

0 

40 

0.32 

1  .  83 

1  .  8976 

500 

1  25 

0 

5 

0.32 

1  .  83 

1.6106 

500 

1  25 

0 

40 

0.32 

3.72 

3.8271 

1000 

125 

0 

5 

0.32 

3.41 

3.3106 

1000 

125 

0 

40 

0.32 

7 . 64 

7 .7848 

250 

250 

0 

5 

0.32 

0 .98 

1 . 0000 

250 

250 

0 

40 

0.32 

0.98 

1  .  0000 

500 

250 

0 

5 

0.32 

1  .  25 

1.2150 

500 

250 

0 

40 

0.32 

1  .87 

1 . 9460 

1000 

250 

0 

5 

0.32 

2.05 

2.0147 

1000 

250 

0 

40 

0.32 

3.82 

3.9110 

250 

500 

0 

5 

0.32 

0.85 

0.9194 

250 

500 

0 

40 

0.32 

0 . 58 

0 . 5376 

500 

500 

0 

5 

0.32 

0.98 

1  .  0000 

500 

500 

0 

40 

0.32 

1  .  00 

1 . 0000 

1000 

500 

0 

5 

0.32 

1 .37 

1 .3540 

1000 

500 

0 

40 

0.32 

1  .  92 

1 .9723 

250 

1000 

0 

5 

0.32 

0.76 

0.8426 

250 

1000 

0 

40 

0.32 

0.36 

0.2917 

500 

1000 

0 

5 

0.32 

0.81 

0.8612 

500 

1000 

0 

40 

0.32 

0.54 

0.5195 
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Two  Dimensional  Model  Output 

Time  Step  =  0.003  seconds 
Convergence  Criteria  =  0.001 


UO.cm/sec  U1 , 

,  cm/sec 

Theta 

1000 

1000 

0 

1000 

1000 

0 

250 

1  25 

30 

250 

125 

30 

0 

30 

500 

125 

30 

0 

30 

1000 

1  25 

30 

250 

250 

30 

250 

250 

30 

500 

250 

30 

500 

250 

30 

0 

30 

1000 

250 

30 

250 

500 

30 

250 

500 

30 

500 

500 

30 

500 

500 

30 

0 

30 

1000 

500 

30 

250 

1000 

30 

250 

1000 

30 

500 

1000 

30 

500 

1000 

30 

1000 

1000 

30 

1000 

1000 

30 

250 

125 

30 

250 

125 

30 

500 

125 

30 

500 

1  25 

30 

1000 

125 

30 

1000 

125 

30 

250 

250 

30 

250 

250 

30 

500 

250 

30 

500 

250 

30 

1000 

250 

30 

1000 

250 

30 

250 

500 

30 

250 

500 

30 

500 

500 

30 

500 

500 

30 

1000 

500 

30 

1000 

500 

30 

250 

1000 

30 

250 

1000 

30 

MODEL 

UNIFIED 

um 

D .  cm 

Eff . 

MODEL  EFF 

5 

0.32 

0 . 98 

1 . 0000 

40 

0.32 

0.99 

1 .0000 

5 

0.32 

0 . 87 

1 . 0405 

40 

0.32 

1  .  28 

1  .6996 

0 

0.32 

ERR 

40 

0.32 

2.63 

3.4431 

0 

0.32 

ERR 

40 

0.32 

5 . 58 

6.9082 

5 

0.32 

0 .89 

0.9918 

40 

0.32 

0.72 

0.8711 

5 

0.32 

0.67 

1  .0796 

40 

0.32 

1  .32 

1  .7273 

0 

0.32 

ERR 

40 

0.32 

2.79 

3.4572 

5 

0.32 

0.85 

0 . 9595 

40 

0.32 

0.43 

0 . 4494 

5 

0.32 

0.76 

0.9842 

40 

0.32 

0.74 

0.8664 

0 

0.32 

ERR 

40 

0.32 

1  .  55 

1  .7308 

5 

0.32 

0.81 

0 . 9299 

40 

0.32 

0.30 

0.2305 

5 

0.32 

0.77 

0 . 9241 

40 

0.32 

0.45 

0.4320 

5 

0.32 

0.70 

0 . 9704 

40 

0.32 

0.81 

0 .8660 

5 

1 

0.79 

1.0130 

40 

1 

0.75 

1  .5025 

5 

1 

0.89 

1 . 0825 

40 

1 

0 . 90 

3  .  1  660 

0 

1 

1 .3858 

40 

1 

4.09 

6 .7622 

5 

1 

0.95 

0 . 9973 

40 

1 

0.63 

0 . 9063 

5 

1 

0.97 

1  .0260 

40 

1 

1  .  05 

1 . 6474 

0 

1 

0.53 

1  .  1649 

40 

1 

2.16 

3.3993 

5 

1 

0.90 

0.9863 

40 

1 

0.48 

0.5927 

5 

1 

0.89 

0.9947 

40 

1 

0.66 

0.8803 

0 

1 

0.85 

1.0516 

40 

1 

1.17 

1.7150 

5 

1 

0.87 

0 . 9748 

40 

1 

0.41 

0.4184 
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Two  Dimensional  Model  Output 

Time  Step  =  0.003  seconds 
Convergence  Criteria  =  0.001 


uo. 

cm/sec 

Ui. cm/sec  Theta 

500 

1000 

30 

500 

1000 

30 

1000 

1000 

30 

1000 

1000 

30 

Uo 

=  wind 

vel oc i ty 

Ui 

=  inlet 

vel oc i t y 

the 

ta  =  sa 

mpling  angle 

d  = 

partic 

le  diameter 

D  = 

inlet 

d i ameter 

MODEL 

UNIFIED 

um  D . 

cm 

Eff . 

MODEL  EFF 

5 

1 

0.87 

0 . 9733 

40 

1 

0.45 

0.4862 

5 

1 

0.86 

0.9896 

40 

1 

0 . 69 

0.8685 
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3D  RESULTS  FOR  SQUARE  INLET  OF  1/2  AREA  =  TO  .5555  CM^2 
DIN  =  DIAMETER  OF  EQUAL  AREA  CIRCLE  AND  IS  USED  TO  CALCULATE 
STOKES  NUMBER 


INLET  AREA:  0.56 


DT 

TRAJ 

BIEM 

UNIFIED 

X 

UO 

UI 

THETA 

DP 

SEC 

CVC 

AREA 

Ai 

Ai 

DIFF 

1000 

125 

0 

5 

0.0003 

0.001 

0.112 

1.6158 

1.8193 

-11% 

1000 

1  25 

0 

5 

0.001 

0 . 05 

0.128 

1  .  8483 

1.8193 

2% 

1000 

125 

0 

5 

0.0003 

0.05 

0.191 

2.7491 

1.8193 

51% 

1000 

125 

0 

40 

0 . 0003 

0.001 

0 . 53 

7 . 6309 

7.2619 

5% 

1000 

125 

0 

40 

0.001 

0 . 001 

0.458 

6 . 5928 

7.2619 

-9% 

250 

1  000 

0 

5 

0.001 

0.01 

1  .733 

0.7799 

0.9500 

-18% 

250 

1000 

0 

5 

0 . 0005 

0 . 05 

1  .  928 

0 . 8677 

0 . 9500 

-9% 

250 

1000 

0 

40 

0.001 

0.01 

0.963 

0.4334 

0.3846 

13% 

250 

1000 

0 

40 

0.0005 

0 .01 

0.944 

0.4248 

0.3846 

10% 

250 

1000 

0 

40 

0 .0003 

0.05 

0.929 

0.4183 

0.3846 

9% 

1000 

125 

30 

5 

0.0003 

0.01 

0.053 

0.7574 

1 .3243 

-43% 

1000 

125 

30 

5 

0 . 0005 

0 . 05 

0.056 

0.8134 

1 .3243 

-39% 

1  000 

125 

30 

40 

0 . 0005 

0 . 05 

0 . 488 

7 .0240 

6 . 6768 

5% 

250 

1000 

30 

5 

0.002 

0 . 05 

1  .  654 

0 . 7446 

0 . 9786 

-24% 

250 

1000 

30 

5 

0.0005 

0 . 05 

1  .  933 

0 . 8698 

0 . 9786 

-1  1% 

250 

1000 

30 

40 

0.001 

0.01 

0 . 58 

0.2611 

0.4622 

-44% 

250 

1000 

30 

40 

0.0005 

0.01 

0.928 

0.4177 

0.4622 

-10% 

250 

1000 

30 

40 

0.0003 

0.05 

0.929 

0.4183 

0.4622 

-10% 

500 

50C 

0 

5 

0.0005 

0.05 

0.57 

1 . 0267 

1 .0000 

3% 

500 

500 

30 

5 

0.0005 

0.05 

0.554 

0.9976 

0.9955 

0% 

500 

250 

0 

5 

0.0005 

0 . 05 

0.314 

1.1321 

1 . 0686 

6% 

500 

250 

30 

5 

0.0005 

0.05 

0.305 

1 . 0972 

1.0219 

7% 

500 

125 

0 

5 

0.0005 

0.05 

0.209 

1 .5023 

1  .  1  936 

26% 

500 

125 

0 

5 

0.0005 

0 . 05 

0.172 

1 . 2408 

1  .  1936 

4% 

500 

125 

30 

5 

0.0005 

0 . 05 

0.164 

1.1801 

1 . 0693 

1  0% 

500 

500 

0 

40 

0. 0005 

0.05 

0 . 573 

1 . 0322 

1 . 0000 

3% 

500 

500 

30 

40 

0 . 0005 

0 . 05 

0.528 

0 . 9506 

0 . 8856 

7% 

500 

250 

0 

40 

0.0005 

0.05 

0 . 544 

1 . 9568 

1.8251 

7% 

500 

250 

30 

40 

0.0005 

0 . 05 

0.496 

1  .7873 

1.6177 

10% 

500 

125 

0 

40 

0.0005 

0 . 05 

0.516 

3.7151 

3 . 4445 

8% 

500 

125 

30 

40 

0.0005 

0.05 

0.439 

3 . 1 646 

3 .0644 

3% 

Uo  =  wind  velocity 

Ui  =  inlet  velocity 

theta  =  sampling  angle 

DP  =  particle  diameter 

DT  =  time  step 

cvc  =  convergence  criteria 

TRAJ  AREA  =  area  Ao  in  freestream 

BIEM  Ai  =  aspiration  efficiency  from  model 

UNIFIED  Ai  =  aspiration  efficiency  fr^m  Hangal  (1990) 
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Two  Dimensional  Particle  Tracking  Program 


C  THIS  PROGRAM  COMBINES  THE  2  DIMENSIONAL  BIEM  SOLUTION  PROGRAM 
C  WRITTEN  BY  LIG6ET  AND  LIU  (1983)  WITH  THE  ANALYTICAL  SOLUTION 
C  TO  THE  PARTICLE  EQUATIONS  OF  MOTION  PROPOSED  BY  ALENIUS 
C  (1989)  TO  PLOT  PARTICLE  TRAJECTORIES  INTO  AN  INLET. 

C 

IMPLICIT  REAL=^8(A-H.O-Z) 

DIMENSION  RLNdOO.  100)  .  RN  (  1  00  .  1  00  )  .  F  (  1  00 )  ,  ID  ( 1  00)  , 
$TITLE(20) .CRITX(2) .CRITY(2) 

COMMON  H (100. 100) .RHS(IOO) .N.X(100).Y(100) .XXNOW. 
$YYNOW.P(100) .PN(100).PN1(5).PN2(5),ICT(100). 

$DPDY(600) ,DPOX(600) .A.L.N1 . N2 . UO , D IN . DP . THETA . U I .XCRIT 
$ . YCR I T . DT , XXNEW . YYNEW . XNOW 1 . YNOW 1 
READ  (7.7)  TITLE 
7  FORMAT  (20A4) 

PRINT  7777. TITLE 
7777  FORMAT  (//1X.20A4) 

READ  (7.*)  KEY 
READ  (7,^)  N 
MP=N 
N1  =  1 

IF  (KEY  .EO.  1)  PRINT  6933. N 
6933  FORMATC  NUMBER  OF  NODES  ='.I3/) 

IF  (KEY  .EQ.  1)  PRINT  6938 

6938  FORMATC . KEY  TO  INPUT  DATA  AND  ID  NUMBERS . 

1'/'  X.Y  ARE  THE  NODAL  COORDINATES  '/'  P.PN  ARE  THE 
2POTENTIAL  AND  THE  NORMAL  DERIVATIVE  AT  THE  NODE') 

IF  (KEY  -EQ.  1)  PRINT  6939 

6939  FORMAT  ('  ID:  IDENTIFIER  FOR  THE  NODAL  UNKNOWNS'/. 

1  '  =2  FOR  P  GIVEN'/, 

2  '  =3  FOR  PN  GIVEN'/) 

IF  (KEY  .EQ.  1)  PRINT  6941 

6941  FORMATC  ICT:  NODAL  CORNER  TYPE  IDENTIFIER'/, 

1'  =1  FOR  PN  GIVEN  APPROACHING  NODE’/, 

2'  =2  FOR  PN  GIVEN  LEAVING  NODE’/, 

3' . ALL  MOVEMENTS  IN  CLOCWISE  SENSE . ' /) 


NODAL  DATA  I/O 
6932  L=N-t-1 

IF  (KEY.EQ. 1 )  PRINT  9953 

9953  FORMATC  PO I  NT '  . 7X . ' X-COOD ’  . 6X ,  ' Y-COOD ’  , 5X , 

1 'POTENTIAL' .3X. 'NORM  DERV ' . 3X , ' ID ' . 8X . ’ ICT ' ) 

DO  6946  1=2, L 

READ  (7,»)  X(I) .Y(I) .P(I) .PN(I) , ID(I) . ICT(I) 

IM1=I-1 

6946  IF  (KEY.EQ. 1)  PRINT  9950 , I M 1 , X ( I ) , Y ( I ) , P ( I )  , PN ( I )  . 
$ID(I) . ICT(I) 

9950  FORMAT(I5.5X.F10.4.2X.F10.4,2X,F10.4,2X,F10.4, 2X, 15. 
15X. 15) 


C 
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C  CLEARING  THE  GLOBAL  MATRIX 

C 

DO  17  1  =  1  .N 
DO  16  J=1.L 

16  H(I . J)=0 . 0 

17  F(I)=0.0 

BOUNDARY  INTEGRATION 

X(1)=X(L) 

X(L+1)=X(2) 

Y(1)=Y(L) 

Y(L+1)=Y(2) 

DO  200  1=2, L 
1 11  =  1 
112=1 
PNJP  =  0 . 0 
P JP  =  0 . 0 
DO  140  J=2.L 
IF(I.EQ.J)  GOTO  130 
IF(I.EQ.J+1)  GOTO  135 
IF(J.NE.L)  GOTO  98 
IF(I .EQ. 2)  GO  TO  135 

9  8  R1=DSQRT((X(J+1  )-X(J))*=*2+(Y(J+1  )-Y(J))**2) 

CO=(X(J-H  )-X(J)  )/R1 
SI=(Y(J+1)-Y(J))/R1 
XIA=(Y(J)-Y(I))*SI+(X(J)-X(I) )*C0 
XIB=(Y(J+1)-Y(I))*SI+(X(J+1)-X(I))*CO 
ETA=DABS((Y(I)-Y(J))*CO-(X(I)-X(J) )*SI) 

100  IFCETA.LT. 1 .E-5)  GO  TO  102 
BTN=DATAN (XIB/ETA) 

ATN=DATAN(XIA/ETA) 

105  ASQ=XIA**2+ETA**2 
ALN=DLOG(ASQ) 

BSQ=XIB**2+ETA’^*2 

106  BLN=DLOG(BSQ) 

107  ONE=BSQ*(BLN-1 . ) -ASQ* (ALN- 1 .) 

TWO=XIB*BLN-XIA»=ALN-2  .*(XIB-XIA)+2  .  *ETA*  (BTN-ATN  ) 
SIGNRN=-(X(J)-X(I))*(Y(J+1)-Y(J)  )  +  (X(  J+1  )-X(J)  )=*(Y(  J)- 

$Y(I)) 

PJ=(XIB^(BTN-ATN)-. 5*ETA*(BLN-ALN) )/(XIB-XIA) 
PJ=DSIGN(PJ .SIGNRN) 

XBXA=DABS(XIB-XIA) 

PNJ=TWO»XIB/(2 .*XBXA)-ONE/(4.*XBXA) 

IF(ICT(J) .NE. 1)  GOTO  321 
PN1 (II1)=PN(J) 

F(I-1 )=F(I-1 )+PNJP*PN(J) 

111=111+1 
PNJP=0 . 0 

321  IF(ICT(J)  .NE.2)  GOTO  432 
F(I-1 )=F(I-1 )+PNJ*PN(J) 

PN2(I I2)=PN(J) 

112=112+1 
PNJ  =  0 . 0 


432  RN(I . J)=PJ+PJP 

RLN(I . J)=PNJ+PNJP 

PJP=(ETA*(BLN-ALN)*0,5-XIA*(BTN-ATN) )/(XIB-XIA) 

PJP=DSIGN(PJP.SIGNRN) 

PNJP=ONE/(4.»=XBXA)-TWO*XIA/(2.*XBXA) 

GO  TO  140 
102  BTN=0.0 
ATN=0.0 
GO  TO  105 
130  ALN=0.0 
BTN=0. 0 
ATN=0.0 

XIB=DSQRT( (X(I+1)-X(I))**2+(Y(I+1 )-Y(I))**2) 

XIA=0 . 0 

ETA=0.0 

ASQ=0.0 

BSQ=XIB**2 

GO  TO  106 

135  XIA=-DSQRT((X(I)-X(J))=^*2+(Y(I)-Y(J))**2) 

BTN  =  0 . 0 

ATN  =  0 . 0 

XIB  =  0 . 0 

BLN=0.0 

BSQ=0-0 

ASQ=XIA**2 

ALN=DLOG(ASQ) 

ETA=0.0 
GOTO  107 

140  CONTINUE 

RN(I .2)=RN(I .2)+PJP 
IF(ICT(2) -NE. 1)  GOTO  544 
F(I-1)=F(I-1 )+PNJP*PN1 (1) 

PNJP=0 . 0 

544  RLN(I , 2)=RLN( I .2)+PNJP 

ALPHA=(  (  (X(I)-X(I-1))*’^2+(Y(I)-Y(I-1  )  >=*=*2  + 
$(X(I)-X(  1  +  1  )  )**2  +  (Y(I)-Y(I  +  1  )  )*=^2-(X(I  +  1  )- 
$  X(I-1)  )**2-(Y(I  +  1)-Y(I-1))=^*2)/(2.*DSQRT(((X(I) 
$-X(I-1))*=^2+(Y(I)-Y(I-1))**2)*(  (X(I)-X(I  +  1 
$))**2+(Y(I)-Y(I  +  1))=^^2)))) 

IF(DABS(ALPHA) .LT. 0.9999)  GOTO  141 
ALPHA=3 . 14159265 
GOTO  170 

141  ALPHA=OACOS(ALPHA) 

IF  (DABS(X(I)-X(I-I)) -GT. 1 .E-5)  GOTO  144 
IF(Y(I) .LT.Y(I-I))  GOTO  142 
IF(X(I+1 ) .LT.X( I ) )ALPHA=6. 2831 8531 -ALPHA 
GOTO  170 

142  IF  (X(I+1) .GT.X(I))ALPHA=6. 28318531-ALPHA 
GOTO  170 

144  A=(Y(I-1)-Y(I))/(X(I-1)-X(I)) 

B=Y(I)-A*X(I) 

YA=A*X(I+1)+B 

IF(Y(I)  .LT.  Y(I-1))  GOTO  150 
IF(X(I)  .LT.  X(I-1))  GOTO  146 
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IF(YA 

.LT. 

Y(I  +  1 

) 

GOTO  1 

70 

146 

IF(YA 

.GT. 

Y(I  +  1 

) 

GOTO  1 

70 

150 

IF(X(I 

)  -GT 

.  X(I 

- 

IF(YA 

.GT. 

Y(I  +  1 

) 

GOTO  1 

70 

152 

IF(YA 

.LT. 

Y(I  +  1 

) 

170 

RN(I  .  I 

)=-ALPHA 

200 

CONTINUE 

13 

FORMAT 

(//) 

ALPHA=6 .28318531 -ALPHA 
ALPHA=6 . 2831 8531 -ALPHA 
))  GOTO  152 

ALPHA=6 . 283  18532-ALPHA 
ALPHA=6 .28318531 -ALPHA 


ASSEMBLING  THE  GLOBAL  EQUATIONS 


DO  1000  1=1 .N 
DO  2000  J=1  .N 
IF(I0(J+1)  .EQ. 
IF(ID(J+1)  .EQ. 
2000  CONTINUE 


2) H(I . J)=-RLN(I+1 , J+1) 

3) H(I . J)=RN(I+1 .J+1) 


COMPUTING  R-H-S  VECTOR 


A=0.0 

DO  2001  J=1  .N 

IF(ICT(J)  .NE.  0)PN(J+1)=0.0 
2001  A=A-P(J+1  )=^RN(I  +  1  .  J+1)+PN(J  +  1)*RLN(I  +  1  .J  +  1) 
H(I .L)=A+F(I) 

1000  CONTINUE 


SOLVING  GLOBAL  EQS  FOR  NODAL  UNKNOWNS 
CALL  EQSOL 

ALLOCATING  EACH  SOLUTION  TO  THE  PROPER  NODAL  UNKNOWN 


JJ1  =  1 

DO  8003  J=2.L 

IF(IDCJ)  -EQ.  2)  GOTO  8005 

P( J)=RHS( JJ1 ) 

JJ1=JJ1+1 
GOTO  8003 

8005  PN( J)=RHS(JJ1 ) 

JJ1=JJ1+1 
8003  CONTINUE 

OUTPUT  BOUNDARY  SOLUTION 


WRITE  (6.13) 

WRITE  (6.30) 

30  FORMAT  (' . COMPUTED  RESULTS . '/) 

do  7654  J=2.1 

write  (6.*)  j- 1 , pn ( j) . p ( j ) 

7654  continue 
C 


o  o  o  o  o  o  o 
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PERFORM  TRAJECTORY  CALCULATION 

INPUT  PARTICLE  AND  AIRSTREAM  VARIABLES  AND  CONV.  PARAMETERS 

DO  4000  K=1 , 2 

IF  (K  -EQ.  1)  THEN 
A=1  . 

ELSE  IF  (K  .EQ.  2)  THEN 
A=-1  . 

ELSE 
ENDIF 
CALL  PART 
CRITX(K)=XCRIT 
CRITY(K)=YCRIT 
4000  CONTINUE 

EFF  =(  ((CRITX(1  )-CRITX(2)  )**2+(CRITY(1  )-CRITY(2)  )**2)*=^.  5 
$*(-1  )’^UO)/(UI*DIN) 

WRITE  (3.*)  UO,UI .DP .DIN . THETA. EFF 

PRINT  *.  EFF 

STOP 

END 

LINEAR  EQUATION  SOLVER  SUBROUTINE 

SUBROUTINE  EQSOL 
IMPLICIT  REAL  ^8(A-H.O-Z) 

DIMENSION  M(100) 

COMMON  A(100.100).X(100).N.X1(100).Y(100) .XXNOW. 

$YYNOW.P( 100) .PN(IOO) .PN1 (5) .PN2(5) . ICT(IOO)  . 

$DPDY(600) .DPDX(600) .A1 .L.N1 . N 2 . UO . D I N . DP , THETA . U I 
$ .XCRIT, YCRIT.DT.XXNEW.YYNEW.XNOW1 . YNOW1 
DO  5  1=1 .N 
M(I)=1 
AMAX=A( I . 1 ) 

DO  2  J=2.N 

IF(DABS(A(I , J) )  .LE.  DABS(AMAX))  GOTO  2 
100  AMAX=A(I . J) 

M(I)=J 

2  CONTINUE 

IFCAMAX  .EQ.  0)  GOTO  98 

3  NN=N+1 

DO  4  J=1 ,NN 

4  A(I . J)=A(I . J)/AMAX 
DO  5  IP=1 .N 

IF(IP  .EO.  I)  GOTO  5 
MMM=M( I ) 

ZMULT=A(IP ,MMM) 

DO  6  J=1 ,NN 

IF(J  .NE.  MMM)  GOTO  9 

8  A(IP,J)=0.0 
GOTO  6 

9  A(IP. J)=A(IP. J)-ZMULT*A(I . J) 

6  CONTINUE 

5  CONTINUE 
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DO  7  1  =  1  .N 
NO=M(I) 

7  X(NO)=A(I .NN) 

RETURN 
98  CONTINUE 
STOP 
END 

SUBROUTINE  FOR  PARTICLE  TRAJ.  CALC.,  2D 

SUBROUTINE  PART 
IMPLICIT  REAL  *8(A-H.O-2) 

DIMENSION  DPDXOL(600) .DPDYOL(600) .MAXI (150) 

COMMON  H(100. 100) .RHS(IOO) .N,X(100),Y(100) .XXNOW. 
$YYNOW.P(100) .PN(IOO) .PN1 (5) .PN2(5) . ICT(IOO) , 

$DPDY(600) .DPDX(600) . A.L.NI . N2 . UO . D IN . DP . THETA . U I 
$ . XCR I T . YCR I T . DT . XXNEW . YYNEW . XNOW 1 , YNOW 1 

VARIABLE  IDENTIFICATION 

UO  =  INITIAL  FREE  STREAM  VELOCITY 
G  =  GRAVITY 

DP  =  PARTICLE  DIAMETER,  CM 

AV  =  AIR  VISCOSITY 

DIN  =  INLET  DIAMETER 

DT  =  TIME  STEP  LENGTH 

PDENS  =  PARTICLE  DENSITY 

ADENS  =  AIR  DENSITY 

VR  =  PARTICLE  RELATIVE  VELOCITY 

RE  =  REYNOLDS  NUMBER 

DPDX,  DPDY  =  AIR  VELOCITY  COMPONENTS 

VXP,  VYP  =  PARTICLE  VELOCITY  COMPONENTS 

XX.  YY  =  PARTICLE  LOCATION  COORDINATES 

THETA  =  NOZZLE  ANGLE 

READ  (8,*)  UO, DP. DIN, DT. THETA. UI .XDIN 
READ  (8.*)  CVC 

PRINT  *.  UO. DP. DIN. DT. THETA. UI .XDIN, CVC 

G=-980. 

dt1=dt 

AV=. 000183 

PDENS=1 . 0 

ADENS=. 001205 

DPL=10000.0*DP 

CC=1  +  ((2.*0.07/OPL)=«(1  .  257+(0.4*EXP(-0.55=*=DPL/0.07)))) 
TAU=  (DP**2=*PDENS=*^CC  )  /  (  1  8*AV) 

THETA  1=-. 01 745 27  7 *THETA 
THETA=.01745277*THETA 
GX=G=^SIN  (THETA) 

GY=G*COS(THETA) 

MANY=599 
NMANY=50 
N2  =  0 

IF  (A  .GT.  0.)  THEN 
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93 


YYN0W=XDIN*DIN*SIN(THETA1 ) 

XXNOW=XDIN*DIN=*=COS(THETA) 

ELSE 

YYNOW=  XD IN*D IN*SIN (THETA1) + ( . 67*XD IN*D IN*COS (THETA) *A) 
XXNOW=XD  IN*D  IN*COS  (THETA)  +  (  .  67=^XD  IN*D  IN^S  I N  (THETA)  *A) 
END  IF 
NCAPF=0 

IF  (A  .GT.  0)  NFCAP=0 

NMISS=0 

CONV=9999 

N1  =  1 


PARTICLE  POSITION  AND  VELOCITY  COMPONENTS 


DO  5000  J=1 .NMANY 
NCAP=0 
N1  =  1 
N2=N2+1 

DO  4300  1=1 .MANY 
CALL  INSOLN 

IF  (I  .EQ.  1)  THEN 

VXNOW=DPDX ( I ) + (TAU*GX) 

VYNOW=DPDY ( I ) + (TAU*6Y) 

ENDIF 

VR=  (  (VXNOW-DPDX  ( I )  )  **2+  (VYNOW-DPDY  (  I ) )  =*=*2 )  **0 . 5 
RE=VR*DP*ADENS/AV 
IF  (RE  .LE.  0.5)  THEN 
QRE=1 

ELSE  IF  ((RE  .GT.  0.5)  .AND.  (RE  . LE .  800.))  THEN 
QRE=1+(0. 15*RE**0.687) 

ELSE 

QRE=0.44*RE/24 

ENDIF 

IF  (I  .EQ.  1)  THEN 
DPDXOL(I)=DPDX(I) 

DPDYOL(I)=DPDY(I) 

ELSE 

DPDXOL(I)=DPDX(I-1) 

DPDYOL(I)=DPDY(I-1) 

ENDIF 

if  ((xxnow  .It.  1.5*din)  .and.  (xxnow  .gt.  -1.5*din) 
$  .and.  (yynow  .It.  1.5*din)  .and. 

$  (yynow  .gt.  -1.5*^din))  then 

dt=dt1/10. 
else 

dt=dt1 
endi  f 

U AX= (DPDX ( I ) -DPDXOL ( I ) ) /DT 
UAY=(DPDY(I)-DPDYOL(I) )/DT 
UBX=DPDX ( I ) - ( (UAX-GX) * (TAU/QRE) ) 
UBY=DPDY(I)-((UAY-GY)*(TAU/QRE)) 
VXNEW=(UAX*DT)+UBX-((UBX-VXNOW)* 

$  EXP(-QRE*DT/TAU)) 
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VYNEW=(UAY*DT)+UBY-( (UBY-VYNOW)^ 

$  EXP(-QRE*DT/TAU)) 

XXNEW=XXNOW+(UAX/2*DT**2)+(UBX*DT)- 
$  (  (UBX-VXNOW)=^(TAU/QRE)*(  1-EXP  (-QRE=^DT/TAU)  )  ) 

YYNEW=YYNOW+(UAY/2*DT**2)+(UBY*DT)- 
$  ( (UBY-VYNOW) * (TAU/QRE) * ( 1 -EXP (-QRE*DT/TAU ) ) ) 

IF  (I  -EQ.  1)  THEN 
XNOW1=XXNOW 
YNOW1=YYNOW 
ENDIF 

WRITE  (4.5002)  XXNOW.YYNOW 
5002  FORMAT  (2F15.10) 

IF  (XXNEW  .LT.  -1 . 1*XDIN/2*DIN)  THEN 
NCAP=2 

PRINT  'PARTICLE  EXIT  REAR  BOUNDARY,  MISS' 

GOTO  4301 

ELSE  IF  ((YYNEW  .GT.  1 . 1 *XD IN*D IN)  .OR.  (YYNEW  .LT. 
$  -1 . 1*XDIN*DIN))  THEN 

NCAP=2 

PRINT  *.  'PARTICLE  EXIT  TOP/BOTTOM.  MISS' 

GOTO  4301 

ELSE  IF  ((XXNEW  .LT.  0.0)  .AND.  (XXNOW  . GT .  0.0) 

$  .AND.  (ABS(YYNEW)  .LT.  DIN/2))  THEN 

NCAP=1 

PRINT  *,  'CAPTURE  AT  FACE' 

GOTO  4301 

ELSE  IF  ((ABS(YYNEW)  .LT.  DIN/2)  .AND. 

$  (ABS(YYNOW)  .GT.  DIN/2)  .AND.  (XXNEW  .LT.  0)) 

$  THEN 

NCAP=2 

PRINT  *.  'CONTACT  TUBE.  MISS' 

GOTO  4301 
ELSE 

NCAP=0 
ENDIF 
N  1  =N 1  + 1 
XXNOW=XXNEW 
YYNOW=YYNEW 
VXNOW=VXNEW 
VYNOW=VYNEW 

4300  CONTINUE 

PRINT  *.  'ERROR,  DID  NOT  EXIT  LOOP' 

PRINT  XNOW1 .YNOW1 .N1 

4301  CONTINUE 
MAXI (J)=N1 

IF  (NFCAP  .GT.  0)  THEN 
GOTO  50 

ELSE  IF  (NCAP  .EQ.  1)  THEN 
NFCAP=1 
YCAP=YNOW1 
XCAP=XNOW1 

YYNOW=YNOW1+( .67*XDIN*OIN*COS(THETA) ) 

XXNOW=XNOW 1 - ( ( YNOW1 -YYNOW) *TAN (THETA) ) 

GOTO  4999 
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ELSE  IF  ((NCAP  .EQ.  2)  .AND.  (YYNOW  .LE.  -OIN/2.))  THEN 
YYNOW=YNOW1+( . 1*DIN) 

XXNOW=XNOW1  +  (  (YYNOW-YNOW1  )=^TAN  (THETA)) 

GOTO  4999 

ELSE  IF  ((NCAP  .EQ.  2)  .AND.  (YYNOW  . 6E .  DIN/2.))  THEN 
YYNOW=YNOW1-( . 1*DIN) 

XXNOW=XNOW 1 - ( (YNOW1 -YYNOW) *TAN (THETA) ) 

GOTO  4999 
ELSE 
ENDIF 

50  IF  (NMISS  .GT.  0)  THEN 

GOTO  10 

ELSE  IF  (NCAP  .EQ.  2)  THEN 
NMISS=1 
YMISS=YNOW1 
XMISS=XNOW1 
GOTO  10 

ELSE  IF  (NCAP  .EO.  1)  THEN 
YYNOW=YNOW1+(A*. 1*DIN) 

XXNOW=XNOW 1 + ( (YYNOW-YNOW1 ) *TAN (THETA) ) 

XDIN=XDIN+( . 1*XDIN) 

GOTO  4999 
ELSE 
ENDIF 

10  IF  (NCAP  .EQ.  1)  THEN 

YCAP=YN0W1 
XCAP=XNOW1 

ELSE  IF  (NCAP  .EQ.  2)  THEN 
YMISS=YNOW1 
XMISS=XNOW1 
ELSE 
ENDIF 

CONV=((XMISS-XCAP)**2+(YMISS-YCAP)**2)**. 5 
IF  (YCAP  .GT.  YMISS)  THEN 

YYNOW=(CONV/2=^COS(THETA)  )+YMISS 
XXNOW=(CONV/2*SIN(THETA))+XMISS 
ELSE  IF  (YMISS  .GT.  YCAP)  THEN 
YYNQW=(CONV/2*COS(THETA))+YCAP 
XXNOW=(CONV/2*SIN(THETA) )+XCAP 
ELSE 
ENDIF 

IF  (CONV  .LT.  CVC)  THEN 
XCRIT=XXNOW 
YCRIT=YYNOW 
GQTO  5001 
ELSE 
ENDIF 

4999  PRINT  *,  N2 . NCAP . YNOW 1 . DPDX ( 1 ) . DPDY ( 1 ) 

5000  CQNTINUE 

5001  PRINT  '  TRAJECTQRIES  HAVE  CONVERGED' 

PRINT  *,  UO,DP,DIN,THETA.XCRIT.YCRIT 
RETURN 

STOP 

END 
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SUBROUTINE  TO  CALC.  VELOCITY  AT  PARTICLE  LOCATION 

SUBROUTINE  INSOLN 
IMPLICIT  REAL  *8(A-H.O-Z) 

DIMENSION  U(600) 

COMMON  H (100, 100) .RHS(IOO) .N.X(100).Y(100) .XXNOW, 
$YYNOW.P(100) .PN(100) .PN1 (5) .PN2(5) , ICT(IOO) . 

$DPDY(600) ,DPDX(600) .A.L.N1 . N2 . UO . D IN . DP . THETA , U I . 
$XCRIT.YCRIT,DT, XXNEW . YYNEW . XNOW 1 . YNOW 1 

PERFORMING  BOUNDARY  INTEG.  FOR  INTERNAL  SOLUTION 

DO  4200  I=N1 .N1 
111  =  1 
112  =  1 
U(I)=0.0 
DPDX(I)=0.0 
DPDY(I)=0.0 
PNJP=0 . 0 
PNXJP  =  0 . 0 
PNYJP=0 . 0 
PJP=0.0 
PXJP=0.0 
PYJP=0.0 
00  4140  J=2.L 

R2=0SQRT((X(J+1)-X(J))**2+(Y(J+1 )-Y(J) )**2) 
CO=(X(J+1)-X(J))/R2 
SI=(Y(J+1 )-Y(J) )/R2 

ETA=DABS( (YYNOW-Y( J))*CO-(XXNOW-X( J) )*SI ) 
XIA=(Y(J)-YYNOW)*SI+(X( J) -XXNOW) *CO 
X IB=  (Y  (J+1)-YYNOW)=^SI  +  (X(J  +  1) -XXNOW)  ^CO 
XBXA=XIB-XIA 

S IGNRN=- (X ( J) -XXNOW) * (Y( J+1 ) -Y ( J) )+(X(J+1 )-X(J) )*(Y(J) 
$  -YYNOW) 

IF(SIGNRN  .LT.  0)ETA=-ETA 
IF(DABS(ETA)  .LT.  1.E-5)  GOTO  4102 
BTN=DATAN(XIB/ETA) 

ATN=DATAN(XIA/ETA) 

GOTO  4105 
4102  BTN=0.0 
ATN=0 . 0 

4105  ASQ=XIA**2+ETA**2 
ALN=DL06(ASQ) 

BSQ=XIB**2+ETA**2 

BLN=DL06(BSQ) 

ONE=BSQ=^(BLN-1  .  ) -ASQ* ( ALN- 1  .  ) 

TWO=XIB*BLN-XIA*ALN-2.*(XIB-XIA)+2  .*ETA=^(BTN-ATN) 
PJ=(XIB*(BTN-ATN)-0.5=^ETA=^(BLN-ALN))/(XIB-XIA) 
PXJ=(-0.5*SI*(BLN-ALN+2.)-CO*(BTN- 
$ATN)+SI*(ETA**2+XIA*XIB)/ASQ)/(XIB-XIA)+ETA*CO/ASQ 
PYJ=(0 . 5*CO*(BLN-ALN  +  2 . )-SI*(BTN-  ATN ) -CO^ (ETA**2+X I  A* 
$XIB)/ASQ)/(XIB-XIA)+ETA*SI/ASQ 
PNJ=TWO*XIB/(2 .*XBXA)-0NE/(4.*XBXA) 
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PNXJ  =  -CO+( (ETA*CO-XIB*SI)*(BTN-ATN)+0 . 5* (X I B*CO+ETA 
$=^SI  )*(BLN-ALN)  )/XBXA 

PNYJ=-SI+( (ETA*SI+XIB*CO)*(BTN-ATN)+0 .5*(XIB*SI- 
$ETA*CO)*(BLN-ALN) )/XBXA 
IF(ICT(J)  .NE.  1)  GOTO  4321 
U(I)=U(I)-PNJP*PN1(I  1 1) 

DPDX(I)=DPDX(I)+PNXJP*PN1 (III) 
DPDY(I)=DPOY(I)+PNYJP*PN1 (II 1) 

111=111+1 
PNJP=0.0 
PNXJP=0 . 0 
PNYJP=0 . 0 

4321  IF(ICT(J)  .NE.  2)  GOTO  4322 
U(I)=U(I)-PNJ*PN2(II2) 

DPDX(I)=DPDX(I)+PNXJ*PN2(1 12) 
DPDY(I)=:DPDY(I)+PNYJ*PN2(II2) 

I.I2  =  I  12+1 

PNJ=0 . 0 
PNXJ=0.0 
PNYJ=0. 0 

4322  U(I)=U(I)+(PJ+FJP)*P(J)-(PNJ+PNJP)*PN(J) 
DPDX(I)=DPDX(I  )  +  (PXJ+PXJP)*P(J)  +  (PNXJ+PNXJP)=^PN(  J) 
DPDY(I)=DPOY(I)+(PYJ+PYJP)*P(J)+(PNYJ+PNYJP)*PN(J) 
PJP=(ETA*(BLN-ALN)*0. 5-XIA*(BTN-ATN) )/(XIB-XIA) 
PXJP=(0. 5*SI*(BLN-ALN-2. )+CO*(BTN-ATN)+SI*(ETA**2+ 

$XIA*XlB)/BSQ)/(XIB-XIA)-ETA*CO/BSQ 
PYJP=(-0. 5*CO*(BLN-ALN-2.  )+SI*(BTN-ATN)-CO*(ETA=*=*2  + 
$XIA=^XIB)/BSQ)/(XIB-XIA)-ETA*SI/BSQ 
PN JP=ONE/ (4 . *XBXA) -TWO*XI A/ ( 2 . *XBXA) 

PNXJP=CO+(  (-ETA*CO+XIA*SI)*(BTN-ATN)-0.5=^(XIA=^CO+ETA 
$SI ) * (BLN-ALN ) ) /XBXA 

PNYJP=SI  +  (-(ETA»=SI+XIA*CO)*(BTN-ATN)-0. 5=^(XIA*SI-ETA 
$  CO)*(BLN-ALN))/XBXA 
4140  CONTINUE 

U(I)=U(I)+PJP*P(2) 

DPDX(I)=DPDX(I)+PXJP*P(2) 

DPDYd  )=DPDY(I  )+PYJP*P(2) 

IF(ICT(2)  .NE.  1)  GOTO  4195 
U(I)=U(I)-PNJP*PN1  (1) 

DP0X(I)=DPDX(I)+PNXJP*PN1 (1) 
DPDY(I)=DPDY(I)+PNYJP*PN1 (1) 

PNJP=0.0 
PNXJP=0.0 
PNYJP=0 . 0 
4195  CONTINUE 

U(I)=U(I)-(PNJP=^PN(2)) 

DPDX(I)=DPDX(I)+PNXJP*PN(2) 

0PDY(I)=DPDY(I)+PNYJP*PN(2) 

U(I)=U(I)/6. 283185308 
DPDX( I )=DPDX( I )/6. 283185308 
DPDY( I )=DPDY( I )/6. 283185308 
4200  CONTINUE 
RETURN 
END 
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Three  Dimensional  Model  Programs: 


PROGRAM  4PTGQ.F 


THIS  IS  A  PROGRAM  TO  CALCULATE  THE  THREE  DIMENSIONAL 
VELOCITY  FIELD  INTO  LOCAL  EXHAUST  OPENINGS  USING  THE 
BIEM  METHOD.  WRITTEN  BY  MIKE  FLYNN  IN  CONJUNCION 
WITH  CASEY  MILLER  IN  MARCH  1988 


IMPLICIT  REAL*8(A-H .0-Z) 

DIMENSION  PHI  (367) .DPHI (3  67) .ALPHA (3 67) . 

1  ID (367) .PREV(3) . AE (69 1 ) . VX ( 1 00) . VY ( 1 00 ) . VZ ( 1 00 ) .R(367)  . 
2H(367 ,367) .U(367) .X(367) . Y(367) .Z(367) ,NE(691 .4) , ICT(367) . 
3XX(100) . YY(IOO) .ZZ(IOO) .ALTDPHI (691  .3) 

COMMON  H . U . MP , X . Y . Z . NE . XX . YY . ZZ . I , J . K . L . A . B . AE 
OPEN  (5,FILE= ' INPUT.DAT* ) 

OPEN  (3.FILE=' INDAT2') 

OPEN  (6,FILE='PRN') 

READ  (5.*)  N 
READ  (5.*)  M 
LL=N+1 
MP=N 

DO  1001  !=1 .N 

READ  (5,=*')  X(I).  Yd).  Z(I).  PHI(I),  DPHI(I).  ID(I).  ICT(I) 
1001  CONTINUE 

READ  (5.*) ((NE(I . II) . 11  =  1 .4)  .  1  =  1  .M) 

PRINT  NE(691,1).  NE(691.2),  NE(691.3).  NE(691.4) 

CLOSE  (3) 


THE  I  LOOP  GENERATES  VALUES  FOR  THE  COEFFICIENT  MATRIX 
H  AND  THE  KNOWN  COLUMN  VECTOR  R.  DIFFERENT  VALUES  FOR  I 
CORRESPOND  TO  THE  BASE  POINT  BEING  POSITIONED  AT  A 
DIFFERENT  NODE.  VALUES  FOR  THE  UNKNOWN  COLUMN  VECTOR 
U  ARE  ALLOCATED  DEPENDING  ON  THE  ID  CODE  AT  THE  START 
OF  THE  LOOP 


DO  10  1=1, N 
R(I)=0.0 
DO  20  J=1 ,N 
H(I  .  J)=0.0 
20  CONTINUE 
10  CONTINUE 

DO  100  1  =  1, N 
IF(ID(I) .EQ. 1)  THEN 
U(I)=DPHI  (I) 

ELSE 

U(I)=PHI  (I) 

ENDIF 

ALPHA(I)=0.0 


ooooooooo  ooooo 
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THE  J  LOOP  CONTAINS  THE  SUMMATIONS  OVER  EACH  ELEMENT 
INITIALLY  THE  AREA  OF  THE  JTH  ELEMENT  IS  CALCULATED 
THIS  IS  AE(J) . 

DO  200  J=1 .M 

DY21=Y(NE('J.2))-Y(NE(J.  1)) 

DY31=Y(NE(J.3))-Y(NE(J,  1)) 

D221=2(NE(J. 2))-Z(NE(J. 1)) 

D231=2(NE(J.3) )-2(NE(J. 1) ) 

DX21=X(NE(J.2))-X(NE(J, 1)) 

DX31=X(NE(J.3))-X(NE(J,  1)) 

AE1=DY21=^D231-D221*DY31 

AE2=D221*DX31-DX21*D231 

AE3=DX21=^DY31-DY21*DX31 

AE(J)=0 . 5*(DSQRT(AE1**2+AE2*=^2+AE3**2) ) 

ASUM=0.0 


THE  K  LOOP  CALLS  THE  PROGRAMS  THAT  PERFORM  THE  BOUNDARY 
INTEGRATIONS  FOR  A  SINGLE  TRIANGULAR  ELEMENT,  SPINT 
IS  USED  IF  THE  BASE  POINT  COINCIDES  WITH  THE  ELEMENT 
NODE  THAT  IS  BEING  CONSIDERED.  GQUAD  IS  THE  STANDARD 
GAUSSIAN  QUADRATURE.  BOTH  PROGRAMS  RETURN  A  AND  B. 


DO  300  K=1 .3 

IF(I .EQ.NE(J.K))  THEN 

CALL  SPINT 

ELSE 

CALL  GQUAD 
ENDIF 

IF(ID(NE(J,K)) -EQ. 1 . AND . ICT (NE ( J , K ) ) .EQ.O)  THEN 
H(I .NE(J.K))=H(I .NE(J.K))-B 
R(I)=R(I)-A*PHI (NE(J.K) ) 

ELSEIF(ICT(NE(J.K)) . EQ . 1 . AND . ID (NE ( J . K ) ) . EQ . 1 . AND . 
$  ID(NE(J.4)) .EQ.2)  THEN 
H(I .NE(J.K))=H(I .NE(J.K)) 

R(  I)=R(I)+B*OPHI  (NE(J,4))-A=^PHI  (NE(J.K)  ) 

ALTDPHI (J.K)=DPHI (NE(J.4) ) 

ELSEIF(ICT(NE(J.K)) .EQ. 1 . AND . ID (NE ( J , K ) ) . EQ . 1 . AND . 
$ID(NE(J.4)) .EQ. 1)  THEN 

H(I .NE(J,K))=H(I .NE(J.K))-B 
R(I)=R(I)-A*PHI (NE(J,K)) 

ELSE 

H(I .NE(J.K))=H(I ,NE(J.K))+A 
R(I)=R(I)+B=^DPHI  (NE(J.K)) 

ENDIF 

ASUM=ASUM+A 

C 

C 

300  CONTINUE 

ALPHA(I)=-A£UM+ALPHA(I) 

200  CONTINUE 


ooo  oooo  a\  oooooo 
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IF(ID(I) .EQ. 1)  THEN 
R  (  I)  =R  (  I  )  -ALPHA  (  I  )  =^PH  1(1) 

ELSE 

H(I . I)=H(I . I)+ALPHA(I) 

ENDIF 

100  CONTINUE 

DO  111  1  =  1  .N 
H(I ,LL)=R(I) 

111  CONTINUE 
CALL  EQSOL 

AFTER  RETURNING  FROM  EQSOL  U(I)  HAS  BEEN  WRITTEN 
OVER  WITH  THE  SOLUTIONS.  THESE  ARE  THEN  ALLOCATED  TO 
PHI (I)  OR  DPHI(I)  AS  FOLLOWS; 

OPEN  (4,FILE='B0UND.0UT‘ ) 

DO  1000  1=1 .N 
IF(ID(I) .EQ. 1)  THEN 
DPHI (I)=U(I) 

ELSE 

PHI (I)=U(I) 

ENDIF 

WRITE  (6.678)  X ( I ) . Y ( I ) . Z ( I ) . PH  I ( I ) . DPH I ( I ) . I D ( I ) , ICT ( I ) 

78  FORMAT  (F 1 5 . 5 .  '  .  ' . F 1 5 . 5 .  '  .  ’  . F 1 5 . 5 ,  '  .  '  . F 1 5 . 5 .  '  .  '  . F  1  5 . 5  ,  ’  .  '  . 

$15.  •  .  M5) 

1000  CONTINUE 

CLOSE  (4) 

READ  IN  THE  INTERNAL  POINTS  FOR  VELOCITY  CALCULATIONS 

READ  (5.*)  NN 
DO  222  1=1 .NN 

READ  (5.*)  XX(I).  YY(I).  ZZ(I) 

222  CONTINUE 

DECLARE  AND  INITIALIZE 

DO  1100  1=1. NN 
PVX=0 . 0 
PVY=0 . 0 
PVZ=0.0 
DO  2000  J=1 ,M 
PREVX=0 . 0 
PREVY=0.0 
PREVZ=0 . 0 
DO  3000  K  =  1  .3 
DO  4000  L  =  1  ,3 
CALL  INTERGQ 

IF(ICT(NE(J.K)) .EQ. 1 .AND. I D (NE ( J  ,  K ) )  .  EQ  .  1 
$.AND. ID(NE(J.4)) .EQ.2)  THEN 
DPHIDE=ALTDPHI (J.K) 

ELSE 

DPHIDE=DPHI (NE(J.K)) 

ENDIF 


oo  o  ooooo 


PREV(L)=(A*PHI (NE(J.K))-B*DPHIDE) 
4000  CONTINUE 

PREVX=PREVX+PREV ( 1 ) 
PREVY=PREVY+PREV(2) 
PREVZ=PREVZ+PREV(3) 

3000  CONTINUE 

PVX=PVX+PREVX 
PVY=PVY+PREVY 
PVZ=PVZ+PREVZ 
2000  CONTINUE 

COEF1=-0. 079577471 
VX(I)=C0EF1*PVX 
VY(I )=C0EF1*PVY 
VZ(I)=C0EF1=^PVZ 
1 1 00  CONTINUE 


OUTPUT  INTERNAL  DATA  TO  FILE 
print  *,  'internal  velocities:' 

OPEN  (4,FILE=' INT.OUT' ) 

DO  444  1=1 .NN 

WRITE  (6,*)  VX(I) , VY(I) .VZ(I) 
444  CONTINUE 

CLOSE  (4) 

CONTINUE 

STOP 

END 


SUBROUTINE  GQUAD 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  W( 13) .CHI (13) ,ETA(13) .F(13) , 
1H(367,367).U(367).X(367).Y(367),Z(367). 

2NE(691 .4) ,XX(100) .YY(IOO) .ZZ(IOO) ,AE(691) 
COMMON  H . U . MP , X . Y , Z , NE . XX . YY , ZZ , I . J . K . L . A , B . AE 
A=0.D00 
B=0 .DOO 

W( 1 )=-0. 56 2 6 00000000000 
W(2)=. 520833333333333 
W(3)=W(2) 

W(4)=W(2) 

CHI  (1 )  =  . 333333333333333 
CHI (2) =.600000000000000 
CHI (3) =. 200000000000000 
CHI  (4)=CHI (3) 

ETA(1)=CHI(1) 

ETA(2)=CHI (3) 

ETA(3)=CHI (2) 

ETA(4)=CHI (3) 

DO  4400  11=1 ,4 
IF(K.EQ.I)  THEN 
F(I I)=CHI (I  I) 

ELSEIF  (K.EQ.2)  THEN 
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F(II)=ETA(II) 

ELSE 

F(II)=1 .D00-ETA(II)-CHI (II) 

ENDIF 

A1  =  (X(NE(J.  1))*CHI(II)+X(NE(J.2))=^ETA(II)+X(NE(J.3)) 
$*(1 -DOO-CHI (II)-ETA(II))-X(I)) 

B1=(Y(NE(J. 1))*CHI (II)+Y(NE(J.2))*£TA(II)+Y(NE(J.3)) 
$*(1 .DOO-CHI (I I)-ETA(I I))-Y(I) ) 

C1=(Z(NE(J. 1))*CHI(II)+Z(NE(J.2))*ETA(I I)+Z(NE(J.3)) 
$*(1 .DOO-CHI (I I)-ETA( II) )-Z(I)) 

R=DSQRT  ( A1  *=*^2+6  1  **2+C  1  **2) 

FF=F(I I)*(1 .DOO/R)*W(I I) 

B=B+FF 

IF(I .EQ.NE(J, 1) .OR. I .EQ.NE(J.2) .OR. I .EQ.NE(J.3))  THEN 
GOTO  4400 
ELSE 

CONTINUE 

ENDIF 

AY1=Y(NE(J. 1))-B1-Y(I) 

AY2=Y(NE( J . 2) )-B1-Y(I) 

AZ1=Z(NE(J, 1))-C1-Z(I) 

AZ2=Z(NE(J,2))-C1-Z(I) 

AX1=X(NE(J, 1))-A1-X(I) 

AX2=X(NE(J,2))-A1-X(I) 

A2=AY1*AZ2-AZ1*AY2 

B2=AZ1=*=AX2-AX1*AZ2 

C2=AX1*AY2-AY1*AX2 

G1=A1*A2+B1*B2+C1*C2 

G2=DSQRT(A2**2+B2**2+C2**2) 

G=F ( 1 1 ) *W ( 1 1 ) * (-G 1 / (R**3*G2 ) ) 

A=A+G 

4400  CONTINUE 
A=AE(J)*A 
B=AE(J)*B 
RETURN 
END 
C 

C  NOTE  THE  SIGN  FOR  G  IS  CORRECT  HERE  ONLY  IF 

C  THE  ELEMENTS  HAVE  BEEN  NUMBERED  COUNTERCLOCKWISE 

C  WHEN  LOOKING  OUT  OF  THE  DOMAIN  THIS  CONVENTION 

C  MUST  BE  OBSERVED  WHEN  CREATING  E(J.K) 

C 

C 

C 

C  THE  SUBRO  TINE  SPINT  IS  DESIGNED  TO  HANDLE  THE  SPECIAL 

C  INTEGRATIONS  THAT  OCCUR  WHEN  AN  ELEMENT  NODE  COINCIDES 

C  WITH  THE  BASE  POINT 

O 

C 

SUBROUTINE  SPINT 
IMPLICIT  REAL*8(A-H.O-Z) 

DIMENSION  W(13) .ETA (13) . CH I ( 1 3) , F( 1 3 ) . 
1H(357.367).U(367).X(367).Y(367),2(367), 

2NE(691 .4)  .XX(IOO)  .YY(IOO)  .ZZdOO)  .AE(691) 


COMMON  H . U . MP . X . Y . Z . NE . XX , YY , 2Z . I . J . K . L . A . B . AE 
A=0 .000 
B=0.D00 
BP  =  0 .000 

W(1) =-0.56 26000000 00000 
W(2)=. 520833333333333 
W(3)=W(2) 

W(4)=W(2) 

CHI (1)=. 333333333333333 
CHI (2) =. 600000000000000 
CHI (3)=. 200000000000000 
CHI (4)=CHI (3) 

ETA(1)=CHI(1) 

ETA(2)=CHI (3) 

ETA(3)=CHI (2) 

ETA(4)=CHI (3) 

DO  500  11=1,4 
IF  (K.EQ.1)  THEN 
F(II)=CHI(II) 

ELSEIF  (K.EQ.2)  THEN 
F(II)=ETA(II) 

ELSE 

F(I I)  =  1  .D00-ETA(I I)-CHI  (I  I) 

ENDIF 

A1=(X(NE(J. 1))*CHI (II)+X(NE(J.2))*ETA(II)+X(NE(J,3)) 
$*(1 .000-CHI (I I)-ETA(ri))-X(I)) 

B1  =  (Y(NE(J,  1  ))=^CHI  (I  I)+Y(NE(J.2)  )*ETA(I  I)+Y(NE(  J.3) ) 
$*(1 .DOO-CHI (II)-ETA(II))-Y(I)) 

C1=(2(NE(J, 1))*CHI (I I)+Z(NE(J.2))*ETA(I I)+2(NE(J.3)) 
$*(1-CHI (II)-ETA(I I))-Z(I)) 

R=OSQRT(A1**2+B1**2+C1**2) 

BPRM31=W(II)*(F(II)-1 .D00)*(1 .DOO/R) 

BP=BP+BPRM31 
500  CONTINUE 

P1P2=DSQRT((X(NE(J.  1))-X(NE(J,2)))=*==*'2  + 

$(Y(NE(J, 1) )-Y(NE(J.2)))**2+(Z(NE(J. 1) )-Z(NE(J. 2)) )**2) 
P1P3=DSQRT(  (X(NE(J.  1))-X(NE(J,3)))=*=*2  + 

$(Y(NE(J, 1))-Y(NE(J.3)))**2+(Z(NE(J. 1))-Z(NE(J,3)))*»2) 
P2P3=DSQRT(  (X(NE(J,2)  )-X(NE(J,3)  ))*»=2  + 

$(Y(NE( J , 2) )-Y(NE(J.3)))**2+(Z(NE(J , 2) )-Z(NE(J,3) ) )**2) 
IF  (K.EQ.1)  THEN 
BS=DABS(P1P3) 

ELSEIF  (K.EQ.2)  THEN 
BS=DABS(P1P2) 

ELSE 

BS=DABS(P2P3) 

ENDIF 

XXX1=0.D00 
YYY2=0 .DOO 
XXX2=BS 

IF(K.EQ.I)  THEN 
BSS=DABS(P1P2) 

ELSEIF(K.EO. 2)  THEN 
BSS=DABS(P2P3) 


o  o  o 
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ELSE 

BSS=DABS(P1P3) 

ENDIF 

IF(K.EQ.I)  THEN 
BSSS=DABS(P2P3) 

ELSEIF(K.EQ.2)  THEN 
BSSS=DABS(P1P3) 

ELSE 

BSSS=DABS(P1P2) 

ENDIF 

THETA1=DACOS( (BS**2+BSS**2-BSSS**2 V(2.D00*BS*BSS) ) 
YYY3=BSS*DS IN (THETA 1) 

XXX3=BSS=^DCOS(THETA1) 

QSQS=DABS(XXX3-XXX2) 

IFCQSQS.LT. 0.0001)  THEN 
PREBS=BSS 
BSS=BS 
BS=PREBS 
XXX2=BS 

YYY3=BSS=^DSIN  (THETA1  ) 

XXX3=BSS*DC0S(THETA1) 

ELSE 

CONTINUE 

ENDIF 

EM=(YYY3-YYY2)/(XXX3-XXX2) 

BB=-EM*BS 

2121  ALP1=(1 .D00/(DSQRT(1 .D00+EM**2))) 

ALP2=  v-EM/ (DSQRT  (  1  .  D00+EM*=^2  )  ) ) 

TNALF2=ALP2/(1 .D00+ALP1) 

TNTHT1=DSIN(THETA1)/(1 .DOO+DCOS(THETA1 )) 
B3NUMR=(TNTHT1+TNALF2)/(1 . DOO-TNTHT 1 *TNALF2 ) 

BSTR3=(  (BB/(DSQRT(1  .D00+EM**2) )  )  =^DL06  (DABS  (B3NUMR/TNALF2 )  )  ) 
B=(AE(J)»‘BP)  +  (BSTR3) 

RETURN 

END 


SUBROUTINE  INTERGQ 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  W(13) . ETA( 1 3) . CH I ( 1 3 ) . F ( 1 3) , 

1H  (3  67 ,367) .U (367) .X(367) . Y (3  67 ) . Z (3  67 )  . 

2NE(691 .4) .XX(IOO) .YY(IOO) .ZZ(IOO) ,AE(691) 
COMMON  H . U . MP . X . Y . Z . NE . XX . YY . ZZ . I . J . K . L . A . B . AE 
A=0.0 
B  =  0.0 

W(1) =-0.562600000000000 
W(2)=. 520833333333333 
W(3)=W(2) 

W(4)=W(2) 

CHI  (1 )  =  .333333333333333 
CHI  (2)  =  .600000000000000 
CHI  (3)  =  . 200000000000000 
CHI  (4)=CHI (3) 


ETA(1)=CHI  (1) 

ETA(2)=CHI (3) 

ETA(3)=CHI (2) 

ETA(4)=CHI (3) 

DO  400  I  1  =  1  .4 
IF(K.EQ.I)  THEN 
F(II)=CHI(II) 

ELSEIF  (K.EQ.2)  THEN 
F(II)=ETA(I I) 

ELSE 

F(II)=1 .-ETA(II)-CHI(II) 

ENDIF 

A1  =  (X(NE(  J,  1)  )=^CHI  (I  I)+X(NE(J.  2)  )*ETA(I  I)+X(NE(J,3) ) 
$*(1-CHI (I I)-ETA(I I) )-XX(I) ) 

B1=(Y(NE(J, 1))*CHI (II)+Y(NE(J.2) )*ETA(I I)+Y(NE(J.3)) 
$=^(1-CHI  (II)-ETA(I  I))-YY(I)) 

C1  =  (Z(NE(J.  1))*CHI  (II)+Z(NE(J.2))=^ETA(I  I)+Z(NE(J.3)) 
$*  (  1  -CH  I  (  I  I )  -ETA  (ID)  -ZZ  (  I) ) 
R=OSQRT(A1**2+B1**2+C1**2) 

IF(L.EO.I)  THEN 
RNUM=A1 

ELSEIF(L.EQ.2)  THEN 
RNUM=B1 
ELSE 
RNUM=C1 
ENDIF 

FF=F(I  I)*(RNUM/R**3)»^W(I  I) 

B=B+FF 

AY1=Y(NE(J. 1))-B1-YY(I) 

AY2=Y(NE(J.2))-B1-YY(I) 

AZ1=Z(NE(J. 1))-C1-ZZ(I) 

AZ2=Z(NE(J. 2) )-C1-ZZ(I) 

AX1=X(NE(J. 1))-A1-XX(I) 

AX2=X(NE(J. 2) )-A1-XX(I) 

A2=AY1*AZ2-AZ1*AY2 
B2=AZ1*AX2-AX1*AZ2 
C2=AX1*AY2-AY1*AX2 
60 1=DSQRT(A2**2+B2**2+C  2=^*2) 

IF(L.EQ-I)  THEN 
CF1=A2 

ELSEIF(L.EQ.2)  THEN 

CF1=B2 

ELSE 

CF1=C2 

ENDIF 

T1=1 . /(R**3*eG1 ) 

T3=T1**2 

T6=A1*A2+B1*B2+C1=*^C2 

6=F  ( I  I )  *W  (  1 1 )  *  (  (-3*RNUM*R*GG  I  *T6*T3 )  +  (CF 1  =^T  1 )  ) 

A=A+G 

400  CONTINUE 

A=AE(J)*A  . 

B=AE(  J)=^B 
RETURN 
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END 

C 

C 

SUBROUTINE  EQSOL 
IMPLICIT  REALMS  (A-H.O-Z) 

DIMENSION  M(3  67) .H (367. 367)  .U(367) 

COMMON  H . U . MP . X . Y , Z . NE . XX . YY . 2Z . I . J . K . L . A . B 
DO  5  11=1 .MP 
M(I I)=1 
AMAX=H(II. 1) 

DO  2  JJ=2.MP 

IF(DABS(H(II . JJ)) -LE.DABS(AMAX))  GOTO  2 
100  AMAX=H ( I  I . JJ) 

M(I I)=JJ 

2  CONTINUE 

IF(AMAX.EQ. 0)  GOTO  98 

3  NN=MP+1 

DO  4  JJ=1 .NN 

4  H(II. JJ)=H(II.JJ)/AMAX 
DO  5  IP=1 .MP 

IF  (IP.EQ. II)  GOTO  5 
MMM=M(I I) 

ZMULT=H(IP.MMM) 

DO  6  JJ=1 ,NN 
IF(JJ.NE.MMM)  GOTO  9 

8  H(IP,JJ)=0.0 
GOTO  6 

9  H(IP. JJ)=H(IP. JJ)-ZMULT*H(II . JJ) 

6  CONTINUE 

5  CONTINUE 

DO  7  11=1 .MP 
NO=M(II) 

7  U (NO)=H ( 1 1 .NN) 

RETURN 

98  PRINT  'no  solution' 

STOP 

END 


o  o  o  o 
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Program  TRK15.F 

C  THIS  PROGRAM  CALCULATES  THE  EFFICIENCY  OF  ASPIRATION 

C  OF  PARTICLES  INTO  AN  ARBITRARY  SHAPED  SAMPLING  INLET 

C  USING  BIEM  TO  CALCULATE  THE  AIR  FLOW  PATTERNS  AND  THE 

C  EQUATIONS  BY  ALENIUS  TO  PREDICT  PARTICLE  LOCATION  AND 

C  AND  VELOCITY  AFTER  A  SMALL  TIME  STEP  HAS  ELAPSED . 

C 

C  THE  PROGRAM  REQUIRES  THE  OUTPUT  OF  THE  BIEM  BOUNDARY 

C  SOLUTION  (4PTGQ.F)  AS  INPUT  IN  ORDER  TO  USE  THE  SHELL 

C  SUBROUTINE  TO  CALCULATE  INTERNAL  VELOCITIES  AT  THE 

C  PARTICLE  LOCATION. 

C 

C  THE  SIZES  OF  THE  ARRAYS  FOR  THE  ARRAY  VARIABLES  IS 

C  DEPENDENT  ON  THE  NUMBER  OF  NODES  AND  ELEMENTS  USED 

C  TO  FORMULATE  THE  BOUNDARY  OF  THE  DOMAIN. 

C 

IMPLICIT  REAL  *8(A-H,0-Z) 

DIMENSION  X(367) .Y(367) . $Z(367) .PHI (367) 

$.DPHI(367) . 10(367) , ICT(367) .NE(693 .4) .AE(693) 

COMMON  UO.UI.DIN.DP.THETAI . DTA , VX . VY . VZ , XX 1 .YY1 ,ZZ1 , 
$NCAP , XX , YY , ZZ , MP , X . Y . Z . NE , I , J . K . L , A . B . AE . XD I N . 

$PHI .DPHI , ID. ICT.N  .M 


INPUT  VARIABLES: 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


uo 

DP 

DIN 

DTA 

THETA1 

UI 

XDIN 

N 

M 

X.Y.Z 
PHI.  DPHI 

ID.  ICT 
NE 


=  WIND  VELOCITY.  CM/SEC 
=  PARICLE  DIAMETER.  CM 
=  INLET  DIAMETER.  CM 
=  INITIAL  TIME  STEP.  SEC 
=  ANGLE  OF  MISALIGNMENT.  DEGREES 
=  INLET  SAMPLING  VELOCITY.  CM/SEC 
=  NUMBER  OF  INLET  DIAMETERS  AWAY  FROM  FACE 
FOR  TRAJECTORY  START  PTS 
=  NUMBER  OF  NODES 
=  NUMBER  OF  ELEMENTS 
=  COORDINATES  OF  NODES 

=  VALUE  OF  POTENTIAL  OR  NORMAL  DERIVATIVE 
0  IF  UNKNOWN 

=  BOUNDARY  CODES.  SEE  4PT6Q.F 

=  ARRAY  NAME  FOR  ALL  THE  TRIANGULAR  ELEMENTS.  EACH 
ELEMENT  IS  IDENTIFIED  BY  THE  NODE  NUMBER  OF 
EACH  CORNER  LISTED  IN  A  CLOCKWISE  DIRECTION 
LOOKING  OUT  OF  THE  DOMAIN.  A  FOURTH  INTERGER 
IN  THIS  ARRAY  SPECIFIES  WHETHER  POTENTIAL  OF 
VELOCITY  IS  KNOWN  OVER  THE  ELEMENT. 


X. Y.Z.PHI .DPHI ,  ID.  AND  ICT  CORRESPOND  TO  THE  OUTPUT  OF 
THE  BIEM  BOUNDARY  SOLUTION  FROM  4PTGQ. 


READ  REQUIRED  INPUT  DATA 


ooo  ro  oo  oooooo 
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READ 
READ 
READ 
READ 
DO  1001 
READ  (5 


(5.*) 

(5.*) 

(5.=^) 

(5.*) 

1  =  1  ,N 
*)  X(I) 


UO.DP.DIN.DTA.THETAI 

CVC.CINC 

N 

M 


UI ,XDIN 


$ICT(I) 

1001  CONTINUE 

READ  (5,*)((NE(I 


Y(I).  Z(I).  PHI(I).  DPHI(I).  ID(I) 


II) . 11=1 .4) , 1=1 .M) 


THETA=THETA1». 017453292 


THE  J  LOOP  CONTAINS  THE  SUMMATIONS  OVER  EACH  ELEMENT 
INITIALLY  THE  AREA  OF  THE  JTH  ELEMENT  IS  CALCULATED 
THIS  IS  AE(J) . 

DO  200  J=1 .M 

DY21=Y(NE(J,2))-Y(NE(J. D) 

DY31=Y(NE(J.3) )-Y(NE(J. 1)) 

D221=Z(NE(J,2))-Z(NE(J,  1)) 

DZ31=Z(NE(J,3) )-Z(NE(J.  1)) 

DX21=X(NE(J.2))-X(NE(J.  1)) 

DX31=X(NE(J.3))-X(NE(J, 1)) 

AE1=DY21*DZ31-DZ21=*DY31 
AE2=D22  1  *0X3  1 -0X2  1  =^DZ3 1 
AE3=DX21*0Y31-DY21*DX31 
AE(J)=0. 5*(OSQRT(AE1**2+AE2**2+AE3**2) ) 

200  CONTINUE 

CALCULATE  FIRST  STARTING  POINT 

NCAP=0 

NTRAJ=0 

XX=XD IN*D IN*COS (THETA1 * (- 1 ) ) 

YY=0.0 

Z2=XDIN*DIN*SIN(THETA1*(-1)) 

50  CONTINUE 

NTRAJ=NTRAJ+1 

IF  (NTRAJ  .EQ.  15)  THEN 

PRINT  *,  "first  Capture  not  found",  NTRAJ 
GOTO  3000 
ENDIF 
CALL  TRAJ 

CHECK  FOR  FIRST  CAPTURE 

IF  (NCAP  .EQ.  1)  THEN 
GOTO  1000 

ELSE  IF  ((NCAP  .EQ.  2)  .AND.  (ZZ  . GT .  0.0))  THEN 
XX=XX1-(CINC*SIN(THETA1) ) 

YY=0.0 

ZZ=ZZ1-(CINC*COS(THETA1)) 

GOTO  250 


ooooo  oooo— too  ooo  ooo 
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ELSE 

XX=XX1  +  (CINC=^SIN(THETA1  )  ) 

YY=0 . 0 

2Z=ZZ1+(CINC*C0S(THETA1 ) ) 

GOTO  250 
ENDIF 

1000  PRINT  "NO.  OF  TRAJ  FCAP  =",NTRAJ 

BEGIN  CRITICAL  TRAJECTORY  CALCULATIONS 

XFCAP=XX1 
YFCAP=YY1 
ZFCAP=ZZ1 

PRINT  ^,"1st  cap xf cap yf cap z f cap 
XCAP=XX1 
YCAP=YY1 
ZCAP=ZZ1 

STEP  =  NUMBER  OF  DEGREES  BETWEEN  RAYS  ON  SEMICIRCLE 
OF  STARTING  POINTS.  NSTEP  =  NUMBER  OF  STEPS 

STEP=15. 

NSTEP=13 
nsum=35 
NCOUNT=0 

INITIAL  CAPTURE  TRAJECTORY  HAS  BEEN  FOUND,  THIS  SECTION  OF 

PROGRAM  CALCULATES  TRAJECTORY  STARTING  POINTS  AT  SUCCESSIVE 
LOCATIONS  AROUND  A  SEMICIRCLE  ABOUT  THE  COORDINATES  OF  THE 
FIRST  CAPTURE 

DO  3000  JJ=1 .NSTEP 
NCOUNT=NCOUNT+1 

theta 2= ( -90+ ((ncount-1)*step))*. 017453292 
AA=.67=^XDIN*DIN 

XX=XFCAP  +  (AA*SIN (THETA2)*S IN (THETA  1  )  ) 

yy=(Aa*cos(theta2) )*(-i; 

ZZ=ZFCAP+(AA*SIN (THETA2)*COS(THETA1 ) ) 

XCAP=XFCAP 
YCAP=YFCAP 
ZCAP=ZFCAP 
NFMISS=0 

ONCE  A  START  POINT  IS  CALCULATED,  THE  KK  LOOP  PLOTS 
SUCCESSIVE  TRAJECTORIES  UNTIL  THE  CONVERGENCE  CRITERIA 
IS  SATISFIED. 

DO  2900  KK=1 ,NSUM 
1500  CALL  TRAJ 

IF  ((NFMISS  .EQ.  0)  .AND.  (NCAP  .EQ.  1))  THEN 
AA=CINC+AA 

XX=XFCAP+(AA*SIN(THETA2)*SIN(THETA1 ) ) 
YY=(AA*C0S(THETA2) )*(-1 ) 


OOOOOO  OOO  OJ  OOOO  tors)  oooo 
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ZZ=ZFCAP+(AA*S1N(THETA2)^C0S(THETA1) ) 
GOTO  1500 

ELSE 

NFMISS=1 

ENDIF 

IF  (NCAP  .EO.  2)  THEN 
XMISS=XX1 
YMISS=YY1 
ZMISS=ZZ1 
ELSE 

XCAP=XX1 

YCAP=YY1 

ZCAP=ZZ1 

ENDIF 


S  IS  THE  DISTANCE  BETWEEN  THE  LAST  CAPTURE  AND  LAST 
MISS.  S  IS  COMPARED  TO  CVC  AND  LOOP  EXITED  IF  S  <  CVC . 


S=SQRT((XCAP-XMISS)**2+(YCAP-YMISS)**2+ 

$  (ZCAP-ZMISS)**2) 

IF  (S  .LT.  CVC)  THEN 
GOTO  2990 
ELSE 

XX=(XMISS+XCAP)/2 

YY=(YMISS+YCAP)/2 

ZZ=(ZMISS+ZCAP)/2 

ENDIF 

if  (kk  .eq.  nsum)  print  *,  "no  convergence,  s  =",s 
900  CONTINUE 

990  XCRIT=(XMISS+XCAP)/2 

YCRIT=(YMISS+YCAP)/2 
ZCRIT=(ZMISS+ZCAP)/2 


EACH  CRITICAL  TRAJECTORY  COORDINATE  IS  WRITTEN  TO  THE 
OUTPUT  UNIT. 


write  (6  "cr  i  t ncount ,  " 

000  CONTINUE 
STOP 
END 


,  xcr i t , " . " , ycri t , " . " , zcr i t 


SUBROUTINE  TRAJ 

THIS  SUBROUTINE  TAKES  A  SINGLE  PARTICLE  STARTING 
POINT.  ASSIGNS  THE  INTIAL  VELOCITY  VALUES  TO  THE 
PARTICLE.  AND  TRACKS  THE  PARTICLE  TO  THE  FACE  OF  THE 
INLET  (CAPTURE)  OR  OUT  OF  THE  BOUNDARY  (MISS) 

IMPLICIT  REAL  *8(A-H.O-Z) 

DIMENSION  VXOL(1600) .VYOL(1600) ,VZOL(1600) . 

IX (367) .Y(367) .Z(367) .PHI (367) .DPHI (367)  .  ID  (367) 

2 . ICT(367) .NE(693 .4) . AE(693) 


1 1 1 


COMMON  UO.UI .DIN.DP.THETA1 . DTA . VX . VY . VZ . XX 1 . YY 1 . ZZ 1 . 

$NC AP , XX . YY . ZZ . MP . X . Y . Z . NE . I . J . K . L , A . B . AE , XD I N . 

$PHI .OPHI . ID. ICT.N.M 
C 

c  print  M 

G=-980 . 

DT=DTA 
AV=. 0001 83 
PDENS=1 . 0 
ADENS=. 001205 
DPL=10000. 0*DP 

CC=1  +  ((2.*0.07/DPL)*(1 .257  +  (0.4=^EXP(-0.5  5=^DPL/0.07)))) 
TAU=(DP**2*PDENS*CC)/(18*AV) 

GX=G*SIN(THETA1 ) 

GZ=G=^COS(THETA1 ) 

GY=0 

C 

C 

MANY=9999 
NCAP=0 
N1  =  1 

DO  4300  1=1 .MANY 
CALL  SHELL 
IF  (I  .EQ.  1)  THEN 
VXNOW=VX+(TAU*GX) 

VYN0W=VY+(TAU*6Y) 

VZNOW=VZ+(TAU*GZ) 

XX1=XX 

YY1=YY 

ZZ1=ZZ 

ENDIF 

VR=( (VXNOW-VX)**2+(VZNOW-VZ)**2 
$  +(VYNOW-VY)**2)**0 . 5 

RE=VR*DP*ADENS/AV 
IF  (RE  -LE.  0.5)  THEN 
QRE=1 

ELSE  IF  ((RE  -GT.  0.5)  .AND.  (RE  . LE .  800.))  THEN 
ORE=1+(0. 15*RE**0 . 687) 

ELSE 

QRE=0 .44*RE/24 
ENDIF 

IF  (I  .EQ.  1)  THEN 
VXOL(I)=VX 
VYOL(I)=VY 
VZOL(I)=VZ 
ELSE 

VXOL(I)=VXOL(I-1) 

VZOL(I)=VZOL(I-1) 

VYOL(I)=VYOL(I-1 ) 

ENDIF 

if  ((XX  .It.  1.5*din)  .and.  (XX  .gt.  (-1.5)=^din) 

$  .and.  (YY  .It.  1.5*din)  .and.  (YY  .gt.  (-1.5)*din) 

$  .AND.  (ZZ  -LT.  1.5=^DIN)  .AND.  (ZZ  .GT.  (-1.5)*DIN)) 

$  then 


DT=DTA/2 . 
el  se 

DT=DTA 
endi  f 

UAX=(VX-VXOL(I))/DT 
UAZ=(VZ-VZOL(I))/DT 
UAY=(VY-VY0L(I))/DT 
UBX=VX-((UAX-GX)*(TAU/QRE) ) 

UBY=VY-( (UAY-6Y)*(TAU/QRE)) 
UBZ=VZ-((UAZ-6Z)*(TAU/QRE) ) 
VXNEW=(UAX*DT)+UBX-((UBX-VXNOW)* 

EXP(-QRE*DT/TAU)) 

VYNEW=CUAY*OT)+UBY-( (UBY-VYNOW)* 

EXP(-QRE*DT/TAU) ) 
VZNEW=(UAZ*DT)+UBZ-((UBZ-VZNOW)* 

EXP(-QRE*DT/TAU)) 

XXNEW=XX+(UAX/2*DT**2)+(UBX*DT)- 

( (UBX-VXNOW) * (TAU/QRE) * ( 1 -EXP ( -QRE*DT/TAU ) ) ) 
ZZNEW=ZZ+(UAZ/2*DT**2)+(UBZ*DT)- 

( (UBZ-VZNOW)=^(TAU/QRE)*(  1-EXP  (-QRE*DT/TAU)  ) ) 
YYNEW=YY+ (UAY/2=*'DT*=^2 )  +  (UBY=^DT) - 

( (UBY-VYN0W)*(TAU/QRE)*(1-EXP(-QRE*DT/TAU) ) ) 

IF  (XXNEW  .LT.  (-1 . 1)*XDIN/2*DIN)  THEN 
NCAP=2 

PRINT  *.  'PARTICLE  EXIT  REAR  BOUNDARY.  MISS' 
GOTO  4301 

ELSE  IF  ((ZZNEW  .GT.  1 . 1 *XD IN*D IN)  .OR.  (ZZNEW  .LT 
(-1  .  D^XDIN’^DIN))  THEN 
NCAP=2 

PRINT  *,  'PARTICLE  EXIT  TOP/BOTTOM.  MISS' 

GOTO  4301 

ELSE  IF  ((YYNEW  . GT .  1 . 1 *XD IN*D I N )  .OR.  (YYNEW  .LT 
(-1.1)  *  XDIN^DIN) )  THEN 
NCAP=2 

PRINT  *.  'PARTRICLE  EXIT  SIDE.  MISS' 

GOTO  4301 

ELSE  IF  ((XXNEW  .LT.  0.05*din)  .AND.  (XX  . GT . 

0.05*din) 

.AND.  (ABS(ZZNEW)  .LT.  DIN/2+.05)  .AND. 

(ABS(YYNEW)  .LT.  DIN/2+.05))  THEN 
NCAP=1 

PRINT  *.  'CAPTURE  AT  FACE' 

GOTO  4301 

ELSE  IF  ((ABS(ZZNEW)  .LT.  DIN/2+.05)  .AND. 
(ABS(YYNEW)  .LT.  DIN/2+.05)  .AND. 

(XXNEW  .LT.  0.0)) 

THEN 

NCAP=2 

PRINT  *.  'CONTACT  TUBE.  MISS' 

GOTO  4301 
ELSE 

NCAP=0 

ENDIF 


oooo  oooooooo  ooo 


XX=XXNEW 

YY=YYNEW 

ZZ=Z2NEW 

VXNOW=VXNEW 

VYNOW=VYNEW 

VZNOW=VZNEW 

4300  CONTINUE 

4301  RETURN 
END 


SUBROUTINE  SHELL 

THIS  SUBROUTINE  CALCULATES  THE  INTERNAL  AIR  VELOCITY 
COMPONENTS  IN  3  DIMENSIONS  AT  THE  LOCATION  THAT  THE  PARTICLE 
IS  LOCATED  USING  THE  BIEM  INTERNAL  SOLUTION  ROUTINE.  IT 
FURTHERC  CALLS  SUBROUTINE  INTERGO  TO  PERFORM  THE  INTERNAL 
GAUSSIAN  QUADRATURE. 


IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  PREV(3) . AE(693) ,NE (693 , 4) , 

IX (367) , Y(367) ,Z(367) , PH  I (367 ) , DPH I (367 )  .  I D (367 ) 

2. ICT(367) .ALTDPHI (693,3) 

COMMON  UO,UI .DIN.OP.THETAl .DTA.VX.VY. VZ.XX1 ,YY1 .ZZ1 . 
$NCAP , XX . YY , ZZ , MP . X , Y . Z , NE . I , J . K . L . A . B . AE . XD I N , 

$PHI .DPHI . ID. ICT.N.M 
MP=N 


DECLARE  AND  INITIALIZE 

PVX=0 . 0 
PVY=0 . 0 
PVZ=0.0 
DO  2000  J  =  1  ,M 
PREVX=0. 0 
PREVY=0 . 0 
PREVZ=0 . 0 
DO  3000  K  =  1  .3 
DO  4000  L=1 ,3 
CALL  INTERGO 

IF(ICT(NE(J.K)) .EQ. 1 .AND. ID(NE(J.K)) .EQ. 1 
$.AND. ID(NE(J.4)) .EQ.2)  THEN 
DPHIDE=DPHI (NE(J.4)) 

ELSE 

DPHIDE=DPHI (NE(J.K)) 

ENDIF 

PREV(L)=(A*PHI (NE(J,K))-B*DPHIDE) 

4000  CONTINUE 

PREVX=PREVX+PREV( 1 ) 

PREVY=PREVY+PREV(2) 

PREVZ=PREVZ+PREV(3) 


o  o  o 
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3000  CONTINUE 

PVX=PVX+PREVX 
PVY=PVY+PREVY 
PVZ=PVZ+PREVZ 
2000  CONTINUE 

COEF1=-0. 079577471 
VX=COEF1*PVX 
VY=COEF1=^PVY 
VZ=COEF1=^PVZ 
1100  CONTINUE 
C 
C 

444  CONTINUE 
C  CLOSE  (4) 

RETURN 

END 


SUBROUTINE  INTERGQ 
IMPLICIT  REAL*8(A-H.O-Z) 

DIMENSION  W(13) ,ETA(13) , CH I ( 1 3) . F ( 1 3) . NE (693 , 4) , 

1X(367) . Y(367) .Z(367) .PHI (367) .DPMI (367) , ID (367) . ICT(367) 
2.AE(693) 

COMMON  UO.UI .DIN.DP.THETAI . DTA . VX . VY , VZ , XX 1 .YY1 .ZZ1 , 

$NC  AP  ,  XX  .  YY  .  ZZ  .  MP  .  X  .  Y  .  Z  .  NE  .  I  .  J  .  K  .  L  .  A  .  B  .  AE  .  XD  I N  . 

$PHI .DPMI . ID. ICT.N.M 
A=0.0 
B  =  0.0 

W(1)=. 333333333333333 
W(2)=W(1) 

W(3)=W(2) 

CHI (1 )=. 6666666666666667 
CHI (2)=. 1 666666666666667 
CHI (3)=CHI (2) 

ETA(1)=. 166666666666667 
ETA(2)=. 166666666666667 
ETA (3) =.666666666666667 
DO  400  11=1,3 
IF(K.EQ.I)  THEN 
F(II)=CHI(II) 

ELSEIF  (K.EQ.2)  THEN 
F(I I)=ETA(II) 

ELSE 

F(II)=1 .-ETA(II)-CHI(II) 

END  IF 

A1=(X(NE(J. 1))*CHI (II)+X(NE(J.2))*ETA(I I)+X(NE(J.3)) 
$*(1-CHI (I I)-ETA(I I))-XX) 

B1=(Y(NE(J. 1 ) )*CHI (II)+Y(NE(J. 2))*ETA(I I)+Y(NE(J.3) ) 
$*(1-CHI (I I)-ETA(II))-YY) 

C1  =  (Z(NE(J.  1 )  )=^CHI  (I  I)+Z(NE(J.  2)  )*ETA(I  I)+Z(NE(J  .3)) 
$*(1-CHI (I I)-ETA(I I))-ZZ) 

R=DSQRT(A1**2+B1**2+C1**2) 

IF(L.EO.I)  THEN 


RNUM=A1 

ELSEIF(L.EO. 2)  THEN 
RNUM=B1 
ELSE 
RNUM=C1 
ENDIF 

FF=F  (ID*  (RNUM/R**3  )  *W  (  I  I  ) 

B=B+FF 

AY1=Y(NE(J, 1))-B1-YY 

AY2=Y(NE(J.2))-B1-YY 

AZ1=2(NE(J. 1))-C1-ZZ 

AZ2=Z(NE(J.2))-C1-2Z 

AX1=X(NE(J. 1))-A1-XX 

AX2=X(NE( J . 2) ) -A1-XX 

A2=AY1*AZ2-AZ1*AY2 

B2=AZ1*AX2-AX1*AZ2 

C2=AX1*AY2-AY1*AX2 

GG1=DSQRT(A2**2+B2**2+C2**2) 

IF(L.EQ.I)  THEN 
CF1=A2 

ELSEIF(L.EQ. 2)  THEN 

CF1=B2 

ELSE 

CF1=C2 

ENDIF 

T1=1 ./(R**3*GG1) 

T3=T1**2 

T6=A1*A2+B1*B2+C1*C2 

G=F(I I)*W(I I)*( (-3*RNUM*R*GG1*T6*T3)  +  (CF1*T1  ) ) 
A=A+G 

400  CONTINUE 
A=AE(J)*A 
B=AE( J) *B 
RETURN 


Program  AREA. FOR 


C  THIS  CODE  CALCULATE  THE  AREA  UNDER  A  CURVE  USING  TRAPIZOIDAL 
RULE. 

C 

C  INPUT 

C  N  ;  NUMBER  OF  INPUT  POINTS 
C  X  ;  X-COORD.  OF  THE  POINT 
C  Y  ;  Y-COORD.  OF  THE  POINT 
C 

C  OUTPUT  - >  AREA  ON  THE  SCREEN. 

C 

JOEBUe 


IMPLICIT  REAL*8(A-H,0-2) 

DIMENSION  X(37) .Y(37) 

READ  (7,*)  N 
DO  100  1  =  1  .N 
READ  (7,*)  y(i),x(i) 

100  CONTINUE 
SUM=0.0 

DO  200  1=1 ,N-1 
0X=X(I+1)-X(I) 

IF(DX.GT.O.OOOOOODO)  THEN 
T=1  .  0 
ELSE 

T  =  -1  .  0 
ENDIF 

TEMP=DABS(DX)*DABS(Y(I+1 )+Y(I) )/2 .DO 
200  SUM=SUM+T*TEMP 

wr 1 te (* , 1 00 1 )  SUM 
1001  format (e20 . 5) 

STOP 

END 


Trajectory  Plotting  Program  in  BASIC 


Program  PART . BAS 


10  OPEN  "filename  of  traj  coords"  FOR  INPUT  AS  #1 
20  OPEN  "flename  of  efficiency  data"  FOR  INPUT  AS  #2 
30  CLS 
40  SCREEN  1 

50  WINDOW  (-6 , -12 )-( 12 , 12  ) 

60  INPUT  #2,  UO,UI, DP, DIN, THETA, EFF 
70  THETA  =  57 . 29577951#*THETA 
80  LINE  (0, DIN/2)-(0, -DIN/2) 

90  LINE  (0,DIN/2)-(-6, DIN/2) 

100  LINE  (0,-DIN/2)-(-6, -DIN/2) 

110  WHILE  NOT  EOF( 1 ) 

120  INPUT  #1 ,  XX, YY 
130  PSET  (XX, YY) 

140  WEND 
150  CLOSE  #1 
160  CLOSE  #2 
170  KEY  OFF 

180  LOCATE  19,8:PRINT  "EFF  =  " : EFF 

190  LOCATE  20, 8: PRINT  "UO  =  " ; UO 

200  LOCATE  2 1,8: PRINT  "DP  =  " ; DP 

210  LOCATE  22,8:PRINT  "DIN  =  ";DIN 

220  LOCATE  23, 8: PR I NT  "THETA  =  "; THETA 

230  LOCATE  24, 8: PRINT  "UI  =  ";UI 

240  A$=INKEY$:IF  A$="C"  THEN  260  ELSE  250 

250  IF  A$="X"  THEN  270  ELSE  240 

260  LPRINT  CHR$(12) 

270  END 
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